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1.  Summary 


We  investigate  the  effects  of  explosion  depth,  tectonic  stresses  and  topography  on  seismic  waves 
from  underground  nuclear  explosions.  We  perform  three-dimensional  calculations  for  an 
explosion  inside  and  under  a  mountain,  and  then  perform  four  three-dimensional  calculations  of 
an  explosion  at  several  depths  in  the  topography  of  the  North  Korean  test  site.  We  also  perform  a 
large  number  of  two-dimensional  axisymmetric  calculations  of  explosions  at  depths  from  150  to 
1000  meters  in  four  earth  structures,  with  compressive  and  tensile  tectonic  stresses  and  with  no 
tectonic  stresses.  We  find  that  P-waves  are  not  strongly  affected  by  any  of  these  effects  because 
the  initial  downgoing  P-wave  is  unaffected  by  interaction  with  the  free  surface.  Surface  waves, 
however,  are  strongly  affected  by  all  of  these  effects.  There  is  an  optimal  depth  where  surface 
waves  are  maximized  at  the  base  of  a  mountain  and  at  or  slightly  below  normal  containment  depth. 
At  deeper  depths,  increasing  overburden  pressure  reduces  the  surface  waves.  At  shallower  depths, 
interaction  with  the  free  surface  reduces  the  surface  waves.  For  explosions  inside  a  mountain, 
displacement  of  the  sides  of  the  mountain  reduces  surface  waves.  Compressive  prestress  reduces 
surface  wave  substantially,  while  tensile  prestress  increases  surface  waves.  The  North  Korean 
explosions  appear  to  be  at  an  optimal  depth,  in  a  region  of  extension,  and  beneath  a  mountain,  all 
of  which  increase  surface  wave  amplitudes.  Measurements  of  Ms  from  Degelen  and  Shagan  River 
explosions  as  small  as  1.5  kilotons  are  consistent  with  the  global  average  Ms:yield  curve,  while 
the  North  Korean  explosions  are  high  by  about  0.6  magnitude  units. 

An  updated  version  of  CRAM3D  is  delivered  along  with  this  report.  The  code  has  an  improved 
procedure  for  bringing  the  grid  into  gravitational  equilibrium,  which  is  particularly  important  for 
calculations  with  topography. 
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2.  Introduction  and  Summary 


The  objective  of  this  project  is  to  decompose  the  seismic  source  to  determine  which  components 
of  the  source  affect  seismic  signals.  We  perform  ID,  2D,  and  3D  simulations  to  identify  the 
physical  mechanisms  operating  during  the  explosion  and  affecting  seismic  signals.  We  focus  the 
data  analysis  on  surface  waves  from  small  (<  50  kt)  events  from  hard  rock  areas  in  Eurasia.  Based 
on  the  modeling  and  data  comparison,  we  assess  the  importance  of  the  components  of  the  explosion 
source  and  their  dependence  on  depth,  material  properties  and  tectonic  environment,  and  determine 
the  implications  for  discrimination  and  yield  estimation. 


We  start  with  data  analysis  by  measuring 
Ms  from  small  explosions  at  the 
Semipalatinsk  test  site.  The  small  event 
analysis  is  motivated  by  the  observation 
that  surface  waves  from  the  North  Korean 
explosions  were  larger  than  expected  for 
their  estimated  yields.  However,  these 
events  were  also  smaller  than  most 
explosions  in  the  global  Ms:yield  data  set 
(e.g.  Stevens  and  Murphy,  2001).  A 
possible  explanation  for  the  anomalous 
yields  is  therefore  that  the  Ms:yield  curve 
has  a  non-unity  slope  for  small  explosions 
which  would  imply  that  Ms  is  generally 
larger  for  smaller  events.  We  test  this 
hypothesis  by  calculating  Ms  for 
explosions  with  yields  as  small  as  1.4 
kilotons.  As  shown  in  Figure  1  and  Figure 


Figure  1.  Data  set  of  new  measurements  together 
with  North  Korean  events. 


2,  Ms  from  these  small  events  are  very  consistent  with  the  global  Ms:  yield  curve,  while  the  North 
Korean  surface  waves  are  well  above  the 


curve. 

We  perform  a  large  number  of  two- 
dimensional  axisymmetric  calculations  of 
explosions  over  a  range  of  depths  from  150 
to  1000  meters  in  four  earth  structures  with 
and  without  tectonic  stresses.  We  calculate 
surface  waves,  body  waves  and  regional 
seismograms  from  all  of  these  calculations. 
We  assess  the  effects  of  depth,  tectonic 
stress  state  and  material  properties  on 
generation  of  seismic  phases.  As  can  be 
seen  in  Figure  3,  the  effect  of  the  nonlinear 
interaction  with  the  free  surface  is  to  reduce 
surface  wave  amplitudes  for  shallow 
explosions.  Surface  waves  also  decrease 
with  increasing  depth  because  of  increasing 
overburden  pressure,  causing  a  peak 
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amplitude  at  about  150  m/kt1/3. 

Compressive  tectonic  stresses  reduce, 
and  sometimes  reverse,  surface  wave 
amplitudes,  while  tensile  tectonic 
stresses  increase  surface  waves, 
although  the  effect  is  smaller.  Also  as 
can  be  seen  by  comparing  Figure  2  with 
Figure  3,  the  new  measurements  from 
Semipalatinsk  are  very  consistent  with 
the  Degelen  calculations,  particularly 
for  a  compressive  stress  state,  while  the 
North  Korean  surface  waves  are  larger 
than  predicted,  even  for  a  tensile 
tectonic  stress  state. 
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Figure  3.  Ms  from  axisymmetric  calculations  in 
Degelen  structure. 


Figure  4  shows  similar  calculations  for 
the  Shoal  material  model  and  earth 
structure.  Because  it  is  a  weaker  granite 

model,  the  surface  waves  are  about  0.2  magnitude  units  larger  for  the  same  yield,  and  at  shallow 
depths  in  tectonic  compression  some  of  the  surface  waves  are  polarity-reversed. 
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We  then  look  at  the  effects  of  topography  by  performing  calculations  of  an  explosion  beneath  a 
sloping  surface  on  a  mountainside. 

This  topic  is  motivated  by  the  fact  that 
many  explosions  occur  in  mountainous 
areas,  and  a  common  way  to  ensure 
containment  is  to  tunnel  into  the  side  of 
a  mountain.  This  has  occurred  at  North 
Korea,  Novaya  Zemlya,  Degelen 
Mountain  and  other  locations.  What 
does  the  presence  of  the  mountain  do  to 
the  seismic  waves  generated  by  the 
explosion  and  to  corresponding  yield 
estimates?  To  address  this  problem,  we 
perform  three  3D  calculations  of  an 
explosion  beneath  a  steep  slope.  The 
first.  Slope  1,  is  placed  beneath  the 
mountain;  the  second,  Slope  2,  is 
higher  and  inside  the  mountain;  while 
the  third,  Slope  3,  is  inside  the 
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Figure  4.  Ms  from  axisymmetric  calculations  in  Shoal 
structure.  Black  circles  indicate  polarity-reversed 
Rayleigh  waves. 


mountain  with  the  yield  increased  by  a  factor  of  4.  All  three  are  performed  using  the  same  earth 
structure  previously  used  for  calculations  of  Shoal  (Stevens  and  Thompson,  2015),  and  the  Shoal 
calculations  provide  a  baseline  reference  model  for  this  series  of  calculations.  We  find  that  surface 
waves  from  explosions  inside  the  mountain  are  substantially  reduced  because  of  horizontal  stress 
release  by  the  sides  of  the  mountain,  however  surface  waves  from  explosions  under  the  mountain 
are  increased. 
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We  perform  four  large  3D  calculations 
of  explosions  at  the  North  Korean  test 
site  using  the  actual  topography  of  the 
test  site  and  at  depths  from  100  meters 
to  800  meters.  Again  we  used  the  Shoal 
material  model.  As  shown  in  Figure  5, 
the  amplification  of  surface  waves 
below  the  mountain  is  substantial,  and 
in  fact  generates  surface  waves  almost 
as  large  as  those  observed  for  North 
Korea.  This  effect  is  primarily  due  to 
the  lower  overburden  pressure  and 
reduced  horizontal  stresses  caused  by 
the  mountain  relative  to  the  same  depth 
in  a  flat  layered  medium.  Comparing 
with  Figure  4,  the  mountain  base 
amplification  effect  with  no  tectonic 
stresses  is  larger  than  the  plane-layered 
medium  with  tensile  tectonic  release, 
and  in  fact  the  physical  mechanism  is 
similar  since  both  produce  reduced 
pressure  and  reduced  horizontal  stresses 


Figure  5.  Ms-log(yield)  vs.  scaled  depth  for  3D  North 
Korea  calculations.  Compare  with  Figure  2.  Multiple 
points  at  each  depth  correspond  to  different  azimuths 
of  observation.  Surface  waves  from  explosions  inside 
the  mountain  are  both  reduced  in  amplitude  and 
increased  in  variability,  while  surface  waves  from 
explosions  beneath  the  mountain  are  amplified. 
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3.  Technical  Approach 
3.1.  The  CRAM3D  Code 


CRAM  3D  is  an  explicit  three-dimensional  Lagrangian  finite  element  code  designed  to  run 
on  multiple  processors  (Stevens  et  al,  2011,  2014a).  For  an  explosion  simulation,  the  cavity  is 
placed  near  the  center  of  the  grid  and  is  enclosed  by  a  spider  grid  which  facilitates  applying  the 
pressure  boundary  condition  and  rezoning  elements  (Figure  6).  The  well-tested  nonlinear  material 
models  from  CRAM  2D  have  been  implemented  in  CRAM  3D.  The  code  includes  gravity  and 
so  includes  the  important  effects  that  result  from  variation  of  overburden  pressure  with 
depth.  Gravitational  equilibrium  is  established  by  running  an  initial  calculation  with  no  source, 
followed  by  a  second  calculation  including  the  explosion  source.  CRAM3D  also  has  the  capability 
to  include  tectonic  prestress  in  the  calculations. 


Figure  6.  The  CRAM  3D  finite  element  outer  grid  (left)  is  rectangular.  The  inner  grid  ( center )  is 
shaped  to  match  the  shape  of  the  explosion  shock  wave.  CRAM2D  uses  a  similar  axisymmetric 
spider  grid  ( right)  in  the  region  around  the  explosion. 


3.2.  Propagation  with  the  Elastodynamic  Representation  Theorem 

The  representation  theorem  allows  us  to  perform  arbitrarily  complex  nonlinear  calculations  in  the 
source  region,  and  then  propagate  them  with  an  appropriate  Green’s  function.  The  representation 
theorem  is  exact.  That  is,  no  matter  how  complex  the  3D  motion  is  on  the  source  region  boundary, 
it  will  be  correctly  propagated  by  the  representation  theorem.  The  only  exception  is  that  it  will  not 
calculate  the  interaction  of  backscattered  waves  reflected  from  outside  the  source  region  with 
complexities  of  the  source  region. 

In  the  three-dimensional  numerical  finite  difference  calculations,  we  save  displacements  and 
stresses  due  to  the  seismic  source  on  a  monitoring  surface  on  the  boundary  of  a  rectangle  (5  planar 
surfaces,  excluding  the  upper  surface),  and  calculate  Green’s  functions  from  each  point  on  the 
monitoring  surface  to  the  receiver  and  so  the  synthetic  seismogram  at  the  receiver  point  X  outside 
of  the  monitoring  surface  is  obtained  by  integrating  over  the  monitoring  surface  sM 

Ui(X)=  j  {(%(£; X)*TJM(0-uy(0*Sik&X)nk}dA 
sM 
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in  the  frequency  domain,  where  G'.(£;  A)  and  S‘jk(^;X)  are  the  Green’s  function  and  the  stress 
tensor  on  the  monitoring  surface  due  to  a  unit  impulsive  force  at  X  in  direction  i,  T™  is  the  traction 

on  the  monitoring  surface  due  to  the  seismic  source,  u  is  the  displacement  on  the  monitoring 
surface,  and  n  is  the  normal  to  the  monitoring  surface.  The  operator  *  denotes  convolution  and  the 
summation  convention  is  assumed. 

We  use  a  plane-layered  Green’s  function  outside  the  source  region.  The  Green’s  functions  for  the 
complete  seismograms  are  derived  from  an  algorithm  based  on  the  work  of  Luco  and  Apsel  (1983) 
and  Apsel  and  Luco  (1983).  The  technique  used  for  surface  waves  is  similar  to  the  method  of 
Bache  et  al.  (1982).  The  Green’s  functions  for  body  waves  are  generated  by  a  procedure  similar  to 
that  described  by  Bache  and  Harkrider  (1976)  using  a  saddle  point  approximation  to  calculate  a 
far-field  plane  wave  for  a  given  takeoff  angle  from  a  source  in  a  plane-layered  medium.  Although 
the  full  waveform  Green’s  functions  generate  the  complete  waveform,  the  other  Green’s  functions 
provide  additional  insight  into  the  source  and  waveform  generation. 

3.3.  Tectonic  Prestress 


5v=Qver 

burden 

Weight 


Normal  ShminOv 
SHmax<Sv 


Performing  a  calculation  of  tectonic  release  requires  an  estimate  of  the  stress  state  in  the  earth. 
While  the  stress  state  varies  between  regions,  considerable  work  has  also  been  done  on  this  subject 
(e.g.,  Zoback,  2010).  In  general  the  vertical  stress  is  equal  to  the  overburden  weight  of  the  material 
above  at  any  given  point.  The  horizontal  stresses  may  be  larger  or  smaller  than  this  value  up  to  the 
point  where  failure  due  to  frictional  sliding  or 
rock  fracturing  relieves  the  stress.  In  most  of 
the  world  tectonic  stresses  are  applied 
continuously,  and  so  the  prestress  is  in  general 
close  to  equilibrium  with  frictional  sliding.  The 
ratio  of  maximum  to  minimum  principal  stress, 
both  reduced  by  pore  pressure,  is  given  by 

Figure  7  Stress  regimes  vary  by  region. 

where  p  is  the  coefficient  of  Horizontal  stresses  may  be  substantially  larger 

or  smaller  than  the  vertical  stress. 
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friction.  For  a  typical  coefficient  of  friction  of 

0.6,  this  ratio  is  3.1,  so  the  difference  between  principal  stresses  can  be  substantial  and  the  effect 
on  seismic  waves  potentially  large.  In  our  calculations,  we  calculate  the  overburden  weight  to 
determine  the  vertical  stress,  and  then  modify  the  horizontal  stresses  to  some  fraction  of  the 
frictional  limit.  This  is  the  initial  stable  state  of  the  calculation  prior  to  introduction  of  the 
explosion.  Note  that  although  the  vertical  stress  is  still  equivalent  to  the  overburden  weight,  the 
pressure  is  not,  and  it  may  be  either  increased  or  reduced  by  the  tectonic  stresses.  Since  material 
strength  increases  with  pressure,  this  also  can  substantially  affect  the  seismic  source.  In  general, 
normal  faulting  regimes  will  amplify  seismic  signals,  while  reverse  faulting  regimes  will  decrease 
seismic  signals;  strike-slip  regimes  could  do  either  (Figure  7). 


3.4  Topography 

CRAM3D  incorporates  topography  by  gradually  distributing  offsets  of  the  grid  from  top  to  bottom 
so  that  change  in  cell  dimensions  is  gradual  rather  than  abrupt  at  the  surface.  The  main  difficulty 
in  implementing  topography  is  establishing  gravitational  equilibrium.  The  variation  in  overburden 
pressure  with  depth  is  very  important  in  nonlinear  calculations  because  strength  is  pressure 
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dependent  and  so  material  becomes  much  weaker  at  shallow  depths.  Because  of  this  we  need  to 
bring  the  grid  into  gravitational  equilibrium  before  running  an  explosion  calculation.  As  mentioned 
above,  this  is  done  by  starting  with  an  approximate  solution  for  overburden  pressure  in  the  grid 
and  then  running  an  initial  calculation  without  the  explosion.  While  this  can  be  accomplished 
relatively  easily  for  a  plane-layered  structure,  bringing  the  grid  to  gravitational  equilibrium  in  a 
structure  with  significant  topography  is  much  more  difficult.  While  the  calculation  will  eventually 
come  to  equilibrium,  adding  topography  is  the  equivalent  of  dropping  a  mountain  onto  the  free 
surface,  and  so  it  can  take  quite  a  long  time  to  come  to  equilibrium. 

The  calculated  displacements  caused  by  topography  close  to  the  NPE  calculation  (Stevens  et  al, 
2014a)  are  shown  in  Figure  8. 


Vertical  Drtlecbon  Caused  by  T errain  {deast=EMmi  Eastward  Deflection  Caused  by  Terrain  (croasl=204m>  Northward  Deflection  Caused  by  Terrain  Cdsast-2ft4m) 


Figure  8.  From  left  to  right,  vertical,  eastward  and  northward  displacements  caused  by 
topography  along  a  slice  204  meters  from  the  NPE  shot  point. 


During  this  project  we  have  made  three  significant  improvements  to  calculations  with  topography 
that  greatly  improve  the  speed  of  the  equilibrium  run  and  also  improve  the  stability  of  equilibrium 
in  the  restart  run.  The  first  improvement  is  to  separate  the  equilibrium  run  into  two  stages.  In  the 
first  stage  we  freeze  the  inner  grid  and  let  the  outer  grid  come  to  approximate  equilibrium  first. 
This  speeds  up  the  equilibrium  run  dramatically  because  the  time  step  for  each  cycle  is  then 
determined  by  the  minimum  grid  size  in  the  outer  grid,  which  is  much  larger  than  the  small  zones 
near  the  cavity  in  the  inner  grid.  In  the  second  stage,  we  release  the  inner  grid  and  let  it  come  to 
equilibrium  with  the  outer  grid.  The  time  step  in  the  second  stage  is  much  shorter,  but  the  outer 
grid  is  already  close  to  equilibrium. 

The  second  improvement  is  to  save  the  body  forces  on  the  bottom  and  sides  of  the  grid  at  the  end 
of  the  equilibrium  run  and  then  to  apply  these  body  forces  in  the  restart  run.  This  guarantees  that 
the  equilibrium  state  will  be  maintained.  Previously  a  more  approximate  solution  was  used  to 
stabilize  the  grid,  but  in  the  project  we  are  modeling  much  larger  topographic  variations,  so  the 
approximate  solution  was  no  longer  adequate. 

The  third  improvement  is  to  use  a  starting  model  that  is  much  closer  to  equilibrium.  See  section 
5.1  for  a  description. 
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3.5  Small  Explosion  Data  Set 

We  acquired  data  from  the  Borovoye  station  for  15  Semipalatinsk  nuclear  tests.  We  used  Russell’s 
(Russell,  2006)  technique  to  measure  Ms.  Table  1  lists  the  parameters  of  these  events  together  with 
the  measured  Ms  at  10  seconds  and  the  maximum  and  minimum  Ms  measured  over  the  recoverable 
frequency  band.  The  best  quality  data  was  near  10s  period. 


T  able  1.  Small  Semipalatinsk  events  with  measurable  Ms  at  Borovoye. 


Event 

Date 

Time 

Yield 

(kt) 

Lat 

Lon 

mb 

Depth 

(m) 

Instrument 

Ms  (10  s) 

Ms(max) 

Ms(min) 

STS_41 4 

12/16/1974 

6:23:00 

3.8 

49.768 

78.082 

4.94 

126 

DS(8-15) 

2.82 

2.94 

2.63 

STS_482 

3/26/1978 

3:57:00 

30 

49.762 

77.983 

5.69 

260 

DS(8-25) 

3.53 

3.61 

3.38 

STS_555 

9/25/1980 

6:21:13 

2.2 

49.783 

78.081 

4.83 

110 

DS(8-13) 

2.51 

2.64 

2.31 

STS_574 

7/17/1981 

2:37:18 

9.3 

49.801 

78.131 

5.07 

146 

DS(8-1 1) 

2.63 

2.75 

2.59 

STS_605 

12/25/1982 

4:23:08 

1.7 

49.781 

78.035 

4.47 

112 

DS(8-25) 

2.24 

2.43 

2.16 

STS_650 

10/18/1984 

4:57:08 

1.4 

49.729 

78.086 

4.25 

106 

DS(9-15) 

2.32 

2.33 

2.23 

STS_676 

5/6/1987 

4:02:08 

32 

49.776 

78.012 

5.6 

174 

SKD(8-25) 

3.68 

3.69 

3.42 

STS_680 

7/17/1987 

1:17:09 

78 

49.776 

78.020 

5.8 

267 

SKD(8-25) 

4.14 

4.18 

4.08 

STS_690 

12/20/1987 

2:55:09 

3.2 

49.776 

78.012 

4.8 

104 

DS(8-13) 

2.48 

2.64 

2.41 

STS_695 

4/22/1988 

9:30:09 

2.3 

49.790 

78.107 

4.9 

124 

DS(10-13) 

2.26 

2.28 

2.18 

STS_698 

6/14/1988 

2:27:09 

5 

50.019 

78.961 

4.8 

271 

DS(8-25) 

2.86 

2.92 

2.83 

STS_701 

9/14/1988 

4:00:00 

140 

49.878 

78.823 

6.03 

651 

DS(8-25) 

4.34 

4.42 

4.30 

STS_703 

11/12/1988 

3:30:06 

17 

50.043 

78.969 

5.24 

0 

DS(8-25) 

3.35 

3.43 

3.29 

STS_704 

11/23/1988 

3:57:09 

19 

49.779 

78.037 

5.4 

204 

DS(9-16) 

3.29 

3.30 

3.24 

STS_71 3 

10/4/1989 

11:29:57 

1.8 

49.748 

78.009 

4.6 

94 

DS(8-16) 

2.58 

2.70 

2.40 
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4.  Results  and  Discussion 


4.1  Analysis  of  Small  Explosions 

As  discussed  in  the  previous  section,  we  measured  Ms  from  Semipalatinsk  explosions  with  yields 
as  low  as  1 .4  kilotons.  The  motivation  for  this  was  to  see  if  the  slope  of  the  Ms:yield  curve  changed 
significantly  for  small  explosions.  The  results  show  that  at  least  for  this  data  set,  it  does  not  (Figure 
9)  and  in  fact  Ms  for  these  events  is  remarkably  consistent  with  the  global  data  set  (Figure  10). 


Figure  9.  Ms  vs.  Yield  for  small  Semipalatinsk  explosions. 


Figure  10.  Ms  vs.  Yield  for  a  global  data  set  of  explosions  from  Stevens  and  Murphy  (2001).  The 
box  shows  the  area  of  the  new  measurements  in  Figure  9. 

In  contrast,  the  Ms:yield  relation  of  the  North  Korean  events  is  very  different.  Figure  1 1  shows  the 
five  North  Korean  nuclear  tests  on  the  same  plot  with  the  small  Semipalatinsk  events.  Here  we 
have  used  the  yield  estimates  from  Murphy  et  al  (2013).  Ms  for  the  North  Korean  events  is  almost 
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a  full  magnitude  larger  at  the  same  yield  compared  to  this  new  data  set.  Ms:mb  shows  a  similar 
anomaly  (Figure  12). 


Figure  11.  Data  set  of  new  measurements  together  with  North  Korean  events. 


Figure  12.  Ms:mb  for  small  Semipalatinsk  and  North  Korean  events. 

Figure  13  shows  mb:yield  for  the  small  Semipalatinsk  events  together  with  the  best  fit  mb:yield 
curve,  which  is  very  close  to  the  mb:yield  curve  from  Murphy  (1995).  There  is  no  mb  anomaly 
comparable  to  the  Ms  anomaly,  although  we  should  state  that  the  yield  estimates  are  derived  from 
mb  and  we  have  no  independent  confirmation  of  yields.  Figure  14  shows  Ms-log(yield)  vs.  depth 
and  scaled  depth  and  there  is  no  apparent  dependence  on  either.  Again,  the  North  Korean  events 
are  very  different. 
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Figure  13.  Left:  mb: Yield  for  small  Semi  events.  Right:  mb: Yield  with  North  Korean  events. 


Depth  (m)  Scaled  Depth  (m/kt1/3) 


Figure  14.  Left:  Ms-log(yield)  vs.  depth.  Right:  Ms-log(yield)  vs.  scaled  depth. 


In  conclusion,  Ms  from  the  small  Semipalatinsk  events  is  very  consistent  with  the  global  data  set 
of  surface  wave  magnitudes  from  explosions,  and  with  an  Ms:log(yield)  slope  of  1.0,  but  very 
inconsistent  with  the  much  larger  surface  waves  of  the  five  North  Korean  explosions. 
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4.2  Numerical  Modeling  of  Seismic  Waves  from  Explosions  Beneath  a  Sloping 
Surface 

Many  explosions  take  place  in  mountainous  areas  and  may  be  adjacent  to  or  beneath  steep  slopes. 
We  model  this  problem  through  a  set  of  3D  finite  difference  calculations.  These  four  calculations 
are  described  in  Table  2.  All  were  performed  in  the  same  uniform  granite  structure  except  for  the 
topographic  deformation.  The  depths  and  yields  were  chosen  to  cover  a  range  of  depths  and  scaled 
depths  that  are  likely  to  be  affected  by  the  sloping  surface. 


Table  2.  Slope  and  Shoal  calculations. 


Name 

Surface 

Depth  Below  Surface 
Vertical/Closest  (m) 

Explosion 
Yield  (kt) 

Scaled  Depth 

m/kt1/3 

Comments 

Shoal 

Flat 

367/  367 

12.5 

158 

Baseline  case 

Slope  1 

Sloped 

458/410 

12.5 

197/177 

2:1  slope 

Slope  2 

Sloped 

207/185 

12.5 

89/80 

2:1  slope 

Slope  3 

Sloped 

207/185 

50 

56/50 

2:1  slope 

The  elevation  profile  used  in  the  slope  calculations  is  shown  in  Figure  15.  Figure  15  also  shows  a 
blow-up  of  the  inner  grid  and  the  final  cavity  within  the  outer  grid  that  is  used  in  the  calculation. 
The  explosion  depth  depends  on  whether  the  depth  is  relative  to  the  closest  point  on  the  surface  or 
the  point  vertically  above  the  explosion.  Both  numbers  are  given  in  Table  2. 


Elevation  (m) 


100 


0 ™0 
0  500  1000  1500  2000  2500  3000 

X  (m) 


Figure  15.  Top:  Elevation  in  map  view  and 
cross  section  for  the  slope  calculation.  Mark 
shows  the  location  of  the  explosion  sources 
in  the  Slope  calculations.  Right:  Inner  and 
outer  grid  near  the  explosion  at  the  end  of  the 
Slope  1  calculation. 


Elevation  (m) 


0  500  1000  1500  2000  2500  3000 


X  (m) 


X  (m) 


The  earth  model  used  with  the  representation  theorem  to  propagate  the  solution  to  regional  and 
teleseismic  distances  is  shown  in  Table  3.  It  is  based  on  a  model  for  the  Climax  stock  area  of  the 
Nevada  Test  Site. 
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Table  3.  Earth  model  used  to  propagate  Shoal  and  Slope  calculations 


Depth 

km 

Thickness 

km 

Vp  km/s 

Vs  km/s 

Density 

g/cm3 

Q 

2.0 

2.0 

5.175 

3.025 

2.600 

100 

3.0 

1.0 

5.500 

3.200 

2.650 

100 

4.0 

1.0 

5.900 

3.400 

2.750 

100 

8.0 

4.0 

5.960 

3.520 

2.780 

350 

12.0 

4.0 

5.960 

3.520 

2.780 

350 

16.0 

4.0 

6.110 

3.610 

2.800 

375 

21.0 

5.0 

6.110 

3.610 

2.800 

375 

31.0 

10.0 

6.370 

3.760 

2.840 

400 

45.0 

14.0 

7.900 

4.420 

3.200 

650 

55.0 

10.0 

8.050 

4.500 

3.300 

600 

65.0 

10.0 

8.050 

4.500 

3.300 

600 

00 

oo 

8.100 

4.500 

3.300 

600 
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4.2.1  Results 

The  figures  on  the  following  pages  show  the  results  of  these  calculations.  They  are  divided  into  4 
sections  -  near  field  results,  regional  full-waveform  seismograms,  fundamental  mode  surface 
waves  and  far-field  body  waves.  The  near  field  results  are  divided  into  two  sections  showing  all 
data  types  for  each  event  and  all  events  for  each  data  type.  Following  is  a  summary  of  results: 

1.  The  reflected  pP  and  pS  waves  are  distorted  relative  to  the  flat  model.  This  causes  a 
reduction  in  the  dominant  pS  phase  although  there  is  additional  S  energy  generated  by  the 
slope.  The  distortion  of  the  pP  phase  causes  a  reduction  in  pP  interference  which  actually 
generates  a  larger  observable  pP  phase. 

2.  Regional  P  and  S  phases  are  slightly  larger  for  the  explosion  in  the  mountain  than  beneath 
it,  especially  for  paths  that  go  across  the  mountain.  Otherwise,  there  is  little  difference  in 
regional  phases  between  the  flat  (Shoal)  and  topographic  (Slope  1  and  Slope  2) 
calculations. 

3.  In  the  Slope  1  calculation  where  the  explosion  is  beneath  the  mountain,  long  period  surface 
waves  show  a  modest  increase  in  amplitude  compared  to  the  flat  (Shoal)  calculation. 

4.  For  Slope  2  and  Slope  3,  where  the  explosion  is  above  the  surrounding  terrain,  surface 
waves  are  reduced  significantly.  This  is  caused  by  a  reduction  in  the  long  period  surface 
wave  due  to  horizontal  stress  relief  by  the  slope.  Stevens  et  al  (1993)  showed  that  this  could 
occur,  although  it  requires  close  proximity  to  a  steep  slope. 

5.  Far  field  SH  waves  which  do  not  exist  for  the  flat  surface,  and  are  quite  small  even  in  the 
presence  of  tectonic  strain  release,  are  comparable  in  amplitude  to  SV  for  the  mountain 
cases. 

Note  that  the  coordinate  system  used  for  the  representation  theorem  calculations  differs  from  the 
coordinate  system  used  for  the  Cram3D  calculations.  In  the  Cram3D  calculations  X  is  east  (along 
the  mountain),  Y  is  north  (across  the  mountain)  and  Z  is  up.  In  the  representation  theorem 
calculations  X  (zero  degrees)  is  north  (across  the  mountain),  Y  (90  degrees)  is  east  (along  the 
mountain)  and  Z  is  down.  See  section  3  of  the  Cram3D  User  Manual  (Stevens  et  al,  2014a)  for  a 
more  detailed  explanation. 


Approved  for  public  release;  distribution  is  unlimited. 

14 


4.2.2  Near  Field  Motion  and  Deformation 


Figure  16.  Shoal  calculation.  Top  row:  vertical  velocity  at  0.1  s  (left),  vertical  velocity  at  0.2  s 
(right).  Bottom  row:  final  plastic  work  (left),  final  crack  strains  (right). 


-1000  -500  0  500  1000  1500  2000  -1000  -500  0  500  1000  1500  2000 

X  (m)  X  (m) 


Figure  17.  Slope  1  calculation.  Top  row:  vertical  velocity  at  0.1  s  (left),  vertical  velocity  at  0.2  s 
(right).  Bottom  row:  final  plastic  work  (left),  final  crack  strains  (right). 
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X  (m)  X  (m) 

Figure  18.  Slope  2  calculation.  Top  row:  vertical  velocity  at  0.1  s  (left),  vertical  velocity  at  0.2  s 
(right).  Bottom  row:  final  plastic  work  (left),  final  crack  strains  (right). 


-1000  -500  0  500  1000  1500  2000  -1000  -500  0  500  1000  1500  2000 


X  (m)  X  (m) 

Figure  19.  Slope  3  calculation.  Top  row:  vertical  velocity  at  0.1  s  (left),  vertical  velocity  at  0.2  s 
(right).  Bottom  row:  final  plastic  work  (left),  final  crack  strains  (right). 
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WORK_slope3.rs,  vz  on:  (y  =  0)  ±  1 
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Figure  20.  Left:  vertical  velocity  at  0.1  seconds  for  the  four  runs.  Right:  vertical  velocity  at  0.2 
second  for  the  four  runs.  From  top  to  bottom:  Shoal,  Slope  1,  Slope  2,  Slope  3. 


1500  2000 
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-1000  -500  0  500  1000  1500  2000  -1000  -500  0  500  1000  1500  2000 

X  (m)  X(m) 

Figure  21.  Left:  Final  plastic  work  for  the  four  runs.  Right:  Final  crack  strains  for  the  four  runs. 
From  top  to  bottom:  Shoal,  Slope  1,  Slope  2,  Slope  3. 
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4.2.3  Regional  Waveforms 

In  this  section  we  show  comparisons  of  regional  waveforms  at  250  km  distance  from  the  Shoal 
and  Slope  calculations.  Figure  22  shows  a  comparison  of  the  regional  waveforms  from  the  Shoal 
and  Slope  1  calculations,  and  a  comparison  with  a  ID  model.  The  ID  calculation  is  a  full  waveform 
seismogram  from  a  point  explosion  source  at  a  depth  of  367  meters  scaled  to  the  moment  of  a  12.5 
kiloton  Mueller-Murphy  (Murphy,  1977)  source  (4.21xl015  N-m).  The  ID  comparison  for  the 
Slope  2  and  Slope  3  calculations  are  for  a  point  source  at  167  meters  depth  (moment  5.48xl015  N- 
m),  and  for  Slope  3  the  moment  is  increased  to  the  moment  of  a  50  kiloton  Mueller-Murphy  source 
(1.83xl016  N-m). 


Figure  22.  Regional  waveforms  at  250  km  from  the  Shoal  calculation  (top)  and  Slope  1 
calculation  (bottom).  Waveforms  from  Slope  1  are  slightly  larger  than  for  Shoal,  and  waveforms 
for  Shoal  are  slightly  smaller  than  the  ID  calculation. 
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Figure  23.  Full  waveforms  at  250  km  from  the  Slope  1  calculation.  Left:  vertical  in  3  directions. 
Right:  Tangential  in  3  directions  and  diagram  of  direction.  (Tangential  at  90  degrees  is  zero). 
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Figure  24.  Full  waveforms  at  250  km  from  the  Slope  2  calculation.  Left:  vertical  in  3  directions. 
Right:  Tangential  in  3  directions  and  diagram  of  direction. 
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Figure  25.  Full  waveforms  at  250  km  from  the  Slope  3  calculation.  Left:  vertical  in  3  directions. 
Right:  Tangential  in  3  directions  and  diagram  of  direction. 
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4.2.4  Surface  waves 


In  this  section  we  show  comparisons  of  long  period  fundamental  mode  surface  waves  (Rayleigh 
and  Love)  at  2000  km  distance  from  the  Shoal  and  Slope  calculations.  Figure  22  shows  a 
comparison  of  Rayleigh  waves  from  the  Shoal  and  Slope  1  calculations,  and  a  comparison  with  a 
ID  model.  The  ID  calculation  is  a  fundamental  mode  Rayleigh  wave  seismogram  from  a  point 
explosion  source  at  a  depth  of  367  meters  convolved  with  a  12.5  kiloton  Mueller-Murphy  source. 
The  ID  comparison  for  the  Slope  2  and  Slope  3  calculations  are  for  a  point  source  at  167  meters 
depth,  and  for  Slope  3  the  yield  is  increased  to  50  kilotons. 
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Figure  26.  Long  period  surface  waves  at  2000  km  from  the  Shoal  calculation  (top)  and  Slope  1 
calculation  (bottom).  Waveforms  from  Slope  1  are  slightly  larger  than  for  Shoal,  and  waveforms 
for  Shoal  are  slightly  smaller  than  the  ID  calculation. 
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Figure  27.  Surface  waves  at  2000  km  for  the  Slope  1  calculation. 
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Figure  28.  Surface  waves  at  2000  km  for  the  Slope  2  calculation. 
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Figure  29.  Surface  waves  at  2000  km  for  the  Slope  3  calculation. 
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4.2.5  Body  waves 


In  this  section  we  show  far  field  P  and  S  body  waves  from  the  Shoal  and  Slope  calculations 
calculated  for  a  uniform  structure  1000  km  below  the  source.  For  comparison,  we  also  show  body 
waves  from  the  Shoal  calculations  with  tensile  and  compressive  prestress  (Stevens  et  al,  2014).  P 
and  SV  waves  are  all  shown  for  a  takeoff  angle  of  30°  while  SH  waves  are  shown  for  three  takeoff 
angles  of  10,  20  and  30  degrees. 
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Figure  30.  Slope  1  P,  together  with  Shoal  calculation  with  no  prestress  and  with  compressive  and 
tensile  prestress.  Similar  plots  follow  for  the  other  cases. 
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Figure  31.  Slope  1  SV. 
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Figure  32.  Slope  1  SH. 
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Figure  33.  Slope  2  P. 
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Figure  34.  Slope  2  SV. 
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Figure  35.  Slope  2  SH. 
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36.  Slope  3  P  (Note  -  this  is  a  50  kt  explosion  compared  to  a  12.5  kt  explosion). 


Velocity  ■  X  Direction  Velocity  -  Y  Direction 


Figure  37.  Slope  3  SV. 
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Figure  38.  Slope  3  SH. 
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4.3  Decomposition  of  the  Seismic  Source  of  Surface  Waves  in  3D  Calculations 

In  the  previous  section  we  looked  at  the  effect  of  source  location  within  and  under  a  mountain. 
Now  we  examine  what  those  variations  do  to  the  seismic  signals.  We  first  look  at  how  surface 
waves  are  generated  by  examining  the  Green’s  functions  and  displacements  on  the  monitoring 
surface.  In  section  3.2,  we  showed  that  the  synthetic  seismogram  at  the  receiver  point  X  outside  of 
the  monitoring  surface  can  be  obtained  by  integrating  over  the  monitoring  surface  sM  using  the 
following  formula: 


Ui(X)=  f  {Gj(|;  X)  *  Tf(0  -  uf  (0  *  Sjfc(£;  X)nk]dA 
Sm 

where  G'(g; X ) and  Sl]k(f\X)  are  the  Green’s  function  and  the  stress  tensor  on  the  monitoring 

surface  due  to  a  unit  impulsive  force  at  X  in  direction  i,  T™  is  the  traction  on  the  monitoring 

surface  due  to  the  seismic  source,  u  is  the  displacement  on  the  monitoring  surface,  and  n  is  the 
normal  to  the  monitoring  surface.  Because  both  the  force  and  impulse  integrated  over  the 
monitoring  surface  must  vanish  in  the  late  time/low  frequency  limits,  the  integral  over  the  Green’s 
function  at  long  periods  is  very  small  (making  sure  this  happens  is  in  fact  one  of  the  more  difficult 
problems  with  numerical  modeling  of  surface  waves  as  a  tiny  amount  of  residual  force  or  impulse 
can  generate  a  large,  spurious  surface  wave).  The  long  period  surface  wave  is  therefore  generated 
by  the  static  displacement  multiplied  by  the  Green’s  function  stresses  integrated  over  the 
monitoring  surface.  Because  of  the  free  surface  boundary  condition,  the  Green’s  function  stresses 
for  vertical  displacements  are  small  at  shallow  depths  and  long  periods,  so  the  surface  waves  are 
primarily  generated  by  the  static  horizontal  displacements.  This  can  be  seen  in  Figure  39  which 
shows  the  numerical  value  of  the  Green’s  function  stresses  at  20  seconds  from  the  surface  to  600 
meters  depth.  The  first  two  indices  in  the  tensor  correspond  to  the  normal  direction  to  the  surface 
and  the  component  of  displacement  at  the  surface  point.  So,  for  example,  in  the  direction  of  the  X 
axis  the  surface  wave  is  primarily  generated  by  the  final  displacement  in  the  X  direction  and 
secondarily  by  the  final  displacement  in  the  Y  direction  with  only  a  small  contribution  from  the 
final  vertical  displacements. 


Figure  39.  Green’s  function  stresses  for  generation  of  a  vertical  component  Rayleigh  wave  (third 
index)  in  the  direction  of  the  X  (left)  and  Y  (right)  axes  at  20  seconds  period.  Long  period 
surface  waves  are  primarily  generated  by  the  static  horizontal  displacements. 
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The  Green’s  functions  on  the  monitoring  surface  are  the  same  for  all  four  of  the  calculations,  but 
the  final  static  displacement  is  not,  and  the  differences  between  them  are  interesting.  Figure  40 
shows  the  horizontal  displacement  normal  to  the  monitoring  surface  in  two  orthogonal  directions. 
The  most  notable  difference  is  the  reduced  displacement  in  Slope  2,  which  leads  to  the  reduced 
surface  wave  amplitudes  shown  earlier.  However,  the  Slope  1  calculation  has  larger  displacement 
than  Shoal,  also  consistent  with  the  surface  waves  shown  earlier.  It  is  not  clear  why  the 
displacement  and  surface  wave  amplitudes  for  Slope  1  are  increased.  One  possible  explanation  is 
that  the  reduced  pressure  under  the  mountain  compared  to  a  flat  surface  at  the  same  depth  causes 
some  additional  dynamic  deformation  near  the  source.  Another  possibility  is  that  the  explosion 
and  mountain  together  cause  more  near-surface  horizontal  displacement  than  a  more  deeply  buried 
source  with  a  flat  surface.  The  effect  of  the  mountain  on  the  near  field  displacements  is  largest 
along  the  direction  of  the  mountain  causing  more  displacement  in  Slope  1  and  less  in  Slope  2  in 
that  direction. 
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Figure  40.  Horizontal  displacements  normal  to  the  monitoring  surface  vs.  depth  for  three 
calculations.  Left:  back  side  of  the  calculation  (across  mountain).  Right:  left  side  of  calculation 
(along  mountain). 
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4.4  Numerical  Modeling  of  Seismic  Waves  from  Explosions  at  the  North 
Korea  Nuclear  Test  Site 

We  perform  calculations  of  nuclear  explosions  at  the  North  Korean  Test  Site,  an  area  with 
substantial  topography  and  where  explosions  have  generated  anomalously  large  surface  waves. 
We  performed  the  four  3-D  calculations  described  in  Table  2  with  the  source  at  a  single  horizontal 
grid  location  but  at  different  depths.  The  best  estimate  of  the  source  depth  of  the  2009  North 
Korean  test  is  540  meters  (Murphy  et  al,  2013),  and  the  other  source  depths  are  selected  to  model 
a  range  of  depths  from  shallower  to  deeper.  An  explosion  yield  of  12.5  kilotons  was  used  in  all 
calculations,  allowing  direct  comparison  with  the  laterally  homogeneous  Shoal  simulation  as  a 
baseline.  Table  2  lists  the  actual  depths  of  burial,  the  minimum  distance  to  surface  and  the  depths 
of  equal  overburden  pressure  in  the  layered  model  bounding  the  grid,  denoted  as  “pressure”  depth 
in  columns  4  and  7.  Note  that  the  pressure  is  lower  in  all  cases  than  the  same  depth  in  a  comparable 
flat  surfaced  model,  and  so  the  “pressure  depth”  is  shallower. 


Table  4.  Summary  of  simulations  of  12.5  kt  explosions  at  North  Korea  test  site. 


Source  Depth  (m) 

Scaled  Source  Depth  (m/kt1/3) 

Name 

Vertical 

Min  to 
Surface 

Pressure 

Vertical 

Min  to 
Surface 

Pressure 

Depth 

Comments 

NK-100 

100 

93 

32 

43 

40 

14 

- 

NK-200 

200 

187 

65 

86 

81 

28 

- 

NK-540 

540 

506 

196 

233 

218 

84 

Estimated  depth 
of  2009  event 

NK-800 

800 

744 

457 

345 

320 

197 

- 

The  grid  was  uniformly  filled  with  the  same  granite  material  model  used  for  the  Shoal  and  Slope 
calculations.  The  regional  earth  model  used  for  propagation  with  the  representation  theorem  is 
specified  in  Table  5.  This  is  the  earth  model  at  the  location  of  the  North  Korean  test  site  from 
Stevens  et  al  (2005),  with  the  2  km  surface  layer  replaced  with  the  properties  of  the  granite  model 
used  in  the  nonlinear  calculations.  Q  was  set  to  a  uniform  value  of  400  in  the  crust  and  upper 
mantle  from  4  to  80  km  depth. 
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Table  5.  Regional  material  model  surrounding  the  North  Korea  test  site 


Depth 

km 

Thickness 

km 

Vp 

km/s 

Vs 

km/s 

Density 

g/cm3 

Q 

2.0 

2.0 

5.175 

3.025 

2.600 

200 

4.0 

2.0 

5.356 

3.100 

2.625 

300 

6.0 

2.0 

5.538 

3.170 

2.650 

400 

8.0 

2.0 

5.719 

3.240 

2.675 

400 

11.0 

3.0 

5.900 

3.312 

2.700 

400 

22.0 

11.0 

6.402 

3.593 

2.736 

400 

32.0 

10.0 

7.016 

3.938 

2.960 

400 

48.0 

16.0 

8.031 

4.508 

3.330 

400 

64.0 

16.0 

7.736 

4.343 

3.223 

400 

80.0 

16.0 

7.431 

4.171 

3.111 

400 

100.0 

20.0 

7.434 

4.173 

3.112 

75 

120.0 

20.0 

7.587 

4.259 

3.168 

75 

142.5 

22.5 

7.733 

4.341 

3.221 

76 

165.0 

22.5 

7.828 

4.394 

3.256 

76 

187.5 

22.5 

7.929 

4.450 

3.293 

77 

210.0 

22.5 

8.110 

4.553 

3.359 

78 

235.0 

25.0 

8.366 

4.696 

3.453 

134 

260.0 

25.0 

8.586 

4.820 

3.533 

135 

00 

OO 

8.742 

4.832 

3.541 

137 

The  simulation  grid  used  the  actual  topography  from  the  North  Korea  test  site.  We  located  the 
source  horizontally  at  the  epicenter  estimated  by  Zhang  and  Wen  (2013)  of  the  2009  test  event, 
illustrated  in  Figure  41,  Panel  (a),  along  with  other  proposed  locations  for  that  event.  The  black 
box  in  Panel  (a)  shows  the  boundaries  of  the  simulation.  This  topography  was  tapered  to  a  uniform 
elevation  to  construct  a  simulation-grid  topography,  shown  in  Panel  (b),  that  allows  for  application 
of  the  representation  theorem  to  propagate  motion  to  farther  distances.  Panels  (c)  and  (d)  show 
south-to-north  (S-N)  and  west-to-east  (W-E)  pressure  profiles  through  the  source  for  an 
isostatically  compensated  earth,  which  was  our  starting  equilibrium.  The  source  is  located  midway 
down  the  long  sloping  south  face  of  the  mountain,  so  that  there  is  significant  topographic 
asymmetry  about  vertical  W-E  plane  through  the  source.  In  contrast  there  is  much  less  asymmetry 
about  the  vertical  S-N  through  the  source.  The  effective  overburden  is  substantially  reduced  within 
or  just  below  the  mountain  compared  to  sources  at  the  same  depth  under  the  layered-earth  structure 
that  bounds  the  grid  as  indicated  in  the  Table  4.  The  “pressure”  depth  under  the  source  epicenter 
is  reduced  by  a  factor  of  2  near  the  surface  and  by  almost  a  factor  of  2  at  800m  actual  depth. 

Grid  nodes  were  spaced  10m  horizontally  ranging  from  1600m  south  to  1400  m  north  and  1700m 
west  to  1700  m  east  of  the  source  epicenter,  and  vertically  every  5m  at  the  grid  edges. 
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Figure  41.  Topography  near  the  North  Korea  test  site,  (a)  shows  the  local  topography  in  the  area 
of  the  test  site  together  with  several  estimated  epicenters  for  the  2009  test.  The  location  of  Zhang 
and  Wen  (2013)  is  used  for  our  simulations.  The  black  line  denotes  the  simulation  area.  Panel 
(b)  shows  the  actual  tapered  topography  used  in  the  simulation.  Panels  (c)  and  (d)  show  south- 
north  and  west-east  pressure  profiles  though  the  source  epicenter. 
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4.4.1  Results 

The  figures  on  the  following  pages  show  the  results  of  these  calculations.  They  are  divided  into  4 
sections  -  near  field  results,  regional  full-waveform  seismograms,  fundamental  mode  surface 
waves  and  far-field  body  waves.  The  near  field  results  include  snapshots  of  velocity  in  N-S  and  E- 
W  cross-section,  and  regions  of  nonlinear  deformation  and  cracking.  Also  included  for  each 
calculation  are  plots  of  outward  final  displacement  on  each  vertical  surface  of  the  monitoring 
surface.  As  noted  in  section  4.3,  these  displacements  are  primarily  responsible  for  surface  wave 
generation. 

Following  is  a  summary  of  results  from  these  calculations: 

1 .  Regional  P  and  S  phases  are  affected  by  the  mountain  topography,  causing  amplification 
relative  to  the  flat-layered  case  in  the  N-S  directions  and  reduction  relative  to  the  flat¬ 
layered  case  in  the  E-W  directions.  For  the  shallowest  two  depths,  regional  P  and  S  waves 
are  reduced  at  azimuths  near  45°,  and  are  as  large  as  or  larger  than  the  flat-layered  case  in 
other  directions.  Average  difference  in  amplitude  is  small,  but  may  be  as  large  as  a  factor 
of  two  for  some  azimuths. 

2.  Long  period  surface  waves  are  significantly  affected  by  both  topography  and  depth.  In 
particular,  surface  waves  are  decreased  substantially  for  shallow  source  depths,  but 
increased  significantly  at  the  base  of  the  mountain.  This  is  strictly  a  long  period  effect, 
easily  visible  when  waveforms  are  low-pass  filtered  at  15  seconds,  but  much  less  apparent 
for  broadband  waveforms. 

3.  Love  waves  are  largest  for  the  100  meter  depth,  and  there  is  also  more  Rayleigh  wave 
variability  for  the  100  meter  depth  case. 

4.  Except  for  the  moveout  in  pP  arrival  time,  far-field  P-waves  are  affected  very  little  by  either 
the  NK  topography  or  source  depth.  There  is  a  slight  reduction  in  P-wave  amplitude  at 
some  azimuths  for  the  shallowest  two  depths. 

5.  Topography  causes  generation  of  far-field  SH  comparable  in  amplitude  to  SV,  and  causes 
more  variability  in  SV.  SV  is  larger  to  the  north  for  the  shallowest  two  depths,  and 
comparable  to  or  smaller  than  SV  from  the  planar  model  for  other  depths  and  directions. 
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4.4.2  Near  Field  Motion  and  Deformation 
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Figure  42.  South-to-north  profiles  (through  the  source)  of  vertical  velocity  from  12.5  kt  explosion 
simulations  at  four  depths  (100m,  200m,  540m,  800m  respectively  from  the  top)  under  sloping 
topography  at  the  North  Korea  test  site  and  the  flat  Shoal  simulation  (bottom):  Left  at  0.1  sec,  right 
at  1.5  seconds. 
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Figure  43.  South- to-north  profiles  (through  the  source)  of  log-magnitude  of  velocity  from  12.5  kt 
explosion  simulations  at  four  depths  (100m,  200m,  540m,  800m  respectively  from  the  top)  under 
sloping  topography  at  the  North  Korea  test  site  and  the  flat  Shoal  simulation  (bottom):  Left  at  0.1 
sec,  right  at  1.5  seconds. 
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Figure  44.  West-to-east  profiles  (through  the  source)  of  vertical  velocity  (left)  and  log-magnitude  of 
velocity  (right)  at  0.2  sec.  from  12.5  kt  explosion  simulations  at  four  depths  (100m,  200m,  540m, 
800m  respectively  from  the  top)  under  sloping  topography  at  the  North  Korea  test  site  and  the  flat 
Shoal  simulation  (bottom). 
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Figure  45.  South-to-north  profiles  (through  the  source)  of  plastic  work  (energy  of  nonlinear 
yielding)  from  12.5  kt  explosion  simulations  at  four  depths  (100m,  200m,  540m,  800m 
respectively  from  the  top)  under  sloping  topography  at  the  North  Korea  test  site,  and  the  Shoal 
simulation  (bottom):  Left  at  0.1  sec,  right  at  1.5  seconds. 

Panel  rows  1-3  on  this  page;  panel  rows  4  and  5  continued  on  the  next  page. 
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Figure  45.  South-to-north  profiles  (through  the  source)  of  plastic  work  (energy  of  nonlinear 
yielding)  from  12.5  kt  explosion  simulations  at  four  depths  (100m,  200m,  540m,  800m 
respectively  from  the  top)  under  sloping  topography  at  the  North  Korea  test  site,  and  the  Shoal 
simulation  (bottom):  Left  at  0.1  sec,  right  at  1.5  seconds.  Panel  rows  4  and  5  on  this  page;  panel 
rows  1-3  on  the  previous  page. 
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Figure  46.  South- to-north  profiles  (through  the  source)  of  crack  strain  from  12.5  kt  explosion 
simulations  at  four  depths  (100m,  200m,  540m,  800m  respectively  from  the  top)  under  sloping 
topography  at  the  North  Korea  test  site  and  the  Shoal  simulation  (bottom):  Left  at  0.1  sec,  right  at 
1.5  sec. 

Panel  rows  1  -3  on  this  page;  panel  rows  4  and  5  continued  on  the  next  page. 
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Figure  46.  South- to-north  profiles  (through  the  source)  of  crack  strain  from  12.5  kt  explosion 
simulations  at  four  depths  (100m,  200m,  540m,  800m  respectively  from  the  top)  under  sloping 
topography  at  the  North  Korea  test  site  and  the  Shoal  simulation  (bottom):  Left  at  0.1  sec,  right  at 
1.5  sec.  Panel  rows  4  and  5  on  this  page;  panel  rows  1-3  on  the  previous  page. 
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Figure  47.  Outward  displacement  on  the  west  monitoring  surface  at  deast  =  -1500m  from  the 
source  for  the  3-D  simulations  of  explosions  at  the  North  Korea  test  site  with  source  depths 
(from  top  to  bottom)  of  100,  200,  540,  and  800  m. 
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Figure  48.  Outward  displacement  on  the  east  monitoring  surface  at  deast  =  1500m  from  the  source 
for  the  3-D  simulations  of  explosions  at  the  North  Korea  test  site  with  source  depths  (from  top  to 
bottom)  of  100,  200,  540,  and  800  m. 
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Figure  49.  Outward  displacement  on  the  north  monitoring  surface  at  deast  =  1500m  from  the 
source  for  the  3-D  simulations  of  explosions  at  the  North  Korea  test  site  with  source  depths 
(from  top  to  bottom)  of  100,  200,  540,  and  800  m. 
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Figure  50.  Outward  displacement  on  the  south  monitoring  surface  at  dnonh  -  -1400m  from  the 
source  for  the  3-D  simulations  of  explosions  at  the  North  Korea  test  site  with  source  depths 
(from  top  to  bottom)  of  100,  200,  540,  and  800  m. 
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Figure  5 1 .  Profiles  of  final  outward  normal  displacement  on  the  monitoring  surface  faces  due 
west,  east,  south,  and  north  of  the  North  Korea  simulation  sources. 
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4.4.3  Regional  Waveforms 

In  this  section  we  show  comparisons  of  regional  waveforms  at  250  km  distance  from  the  North 
Korea  and  Shoal  calculations.  Figure  52  to  Figure  56  show  a  comparison  of  the  regional  waveforms 
from  the  Shoal  and  North  Korea  calculations.  Note  that  both  the  Shoal  and  North  Korea 
calculations  shown  here  were  propagated  with  the  North  Korea  structure  listed  in  Table  5  rather 
than  the  Western  US  structure  used  earlier  for  the  Shoal  and  Slope  calculations.  In  all  of  these 
calculations,  the  azimuth  shown  is  true  azimuth  measured  clockwise  from  North  (Figure  41). 


Approved  for  public  release;  distribution  is  unlimited. 

48 


Regional  (250  km)  V.  :  Zsrc=100m 


Regional  (250  km)  V.  :  Zsrc=200m 


(/) 

E 


> 


l/l 

E 


> 


1/1 

E 


1/1 

E 

> 


Regional  (250  km)  V r  :  Zsrc=800m 


az= 

=0° 

L 

- Shoal 

Hr- 

x  10'4, 

- N.  Korea 

i 

1 

az= 

in 

.. 

i 

JL 

“FT 

- Shoal 

r 

xlO'4 

| 

- N. Korea 

1 

i 

az= 

=-90° 

1 

_ 

i 

l 

- Shoal 

Hr 

*  io 4 

i 

1  - N.  Korea 

az= 

=-135° 

...  .  J 

iL. 

x  10'4 

1 

- Shoal 

- N.  Korea 

ir 

az= 

1 

=-180° 

L 

_ 

- Shoal 

- N.  Korea 

If 

40  50  60  70  80  90  100 

time  (s) 


Figure  52.  Full-waveform  radial  velocity  at  250  km  and  azimuths  of  0,  -45,  -90,  -135  and  -180 
degrees  from  the  North  Korea  simulations  with  source  depths  of  100,  200,  540  and  800m  (red) 
and  the  Shoal  simulations  at  367  m  below  a  flat  surface  (blue). 


Approved  for  public  release;  distribution  is  unlimited. 

49 


Regional  (250  km)  Vf  :  ZJ#r=100m  Regional  (250  km)  V. :  Zjre=200m 


Figure  53.  Full-waveform  vertical  velocity  at  250  km  and  azimuths  of  0,  -45,  -90  -135  and  -180 
degrees  from  the  North  Korea  simulations  with  source  depths  of  100,  200,  540  and  800m  (red) 
and  the  Shoal  simulations  at  367  m  below  a  flat  surface  (blue).  NK  travel  times  were  reduced  by 
amounts  indicated  at  the  top  right  of  each  frame  to  attain  time  alignment  with  the  Shoal 
waveforms.  This  figure  shows  the  early  part  of  the  Figure  52  waveforms. 
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Figure  54.  Full-waveform  tangential  velocity  at  250  km  and  azimuths  of  0,  -45,  -90  -135  and 
-180  degrees  from  the  North  Korea  simulations  with  source  depths  of  100,  200,  540  and  800m 
(red)  and  the  Shoal  simulations  at  367  m  below  a  flat  surface  (blue). 


Approved  for  public  release;  distribution  is  unlimited. 

51 


az=0° 

r 

- Shoal 

f 

xlO'4 

- N.  Korea 

! 

az= 

=-45° 

It 

- Shoal 

r 

x  10"4 

- N.  Korea 

f 

az= 

=-90° 

i 

L 

x  10"4 

- Shoal 

- N.  Korea 

t 

az= 

=-135°  ' 

1 

i 

xlO'4 

- Shoal  ’ 

- N.  Korea 

az= 

, .  i 

=-180° 

L . 

—  Shoal 

r- 

- N.  Korea 

1 _ i _ L 

40  50 


60  70 

time  (s) 


90  100 


Regional  (250  km)  V,  :  Zsrc=540m 


Regional  (250  km)  Vz  :  Zsrc=800m 


Figure  55.  Full-waveform  vertical  velocity  at  250  km  and  azimuths  of  0,  -45,  -90  -135  and  -180 
degrees  from  the  North  Korea  simulations  with  source  depths  of  100,  200,  540  and  800m  (red) 
and  the  Shoal  simulations  at  367  m  below  a  flat  surface  (blue). 
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Figure  56.  Full-waveform  vertical  velocity  at  250  km  and  azimuths  of  0,  -45,  -90  -135  and  -180 
degrees  from  the  North  Korea  simulations  with  source  depths  of  100,  200,  540  and  800m  (red) 
and  the  Shoal  simulations  at  367  m  below  a  flat  surface  (blue).  NK  travel  times  were  reduced  by 
amounts  indicated  at  the  top  right  of  each  frame  to  attain  time  alignment  with  the  Shoal 
waveforms.  This  figure  shows  the  early  part  of  the  Figure  55  waveforms. 
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4.4.4  Surface  Waves 


In  this  section  we  show  comparisons  of  long  period  fundamental  mode  surface  waves  (Rayleigh 
and  Love)  at  2000  km  distance  from  the  Shoal  and  North  Korea  calculations,  again  all  propagated 
with  the  North  Korea  structure  listed  in  Table  5.  Note  that  there  is  only  a  single  Shoal  calculation 
at  370  meters  depth,  so  it  does  not  change  with  the  depth  of  the  North  Korean  calculations.  Figure 
58,  Figure  59  and  Figure  60  show  the  broadband  fundamental  mode  surface  wave  for  the  radial, 
tangential  and  vertical  components,  respectively.  Figure  61,  Figure  62  and  Figure  63  show  the 
same  waveforms  low-pass  filtered  at  15  seconds.  As  with  the  regional  waveforms,  the  azimuth 
shown  is  true  azimuth  measured  clockwise  from  North  (Figure  41).  Figure  57  shows  Ms  calculated 
from  the  vertical  component  for  all  surface  waves.  At  the  shallowest  depth,  there  is  considerable 
variability  with  azimuth.  Surface  waves  from  deeper  events  show  amplification  with  a  peak  at  the 
540  meter  depth.  This  is  most  likely  because  the  pressure  and  horizontal  stress  at  a  depth  beneath 
the  top  of  a  mountain  are  lower  than  at  the  same  depth  in  a  flat  structure. 


Figure  57.  Left:  Ms  vs  depth  for  calculated  surface  waves.  Right:  Ms-log(yield)  vs.  scaled  depth. 
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Figure  58.  Surface  wave  radial  velocity  at  2000km. 
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Figure  59.  Surface  wave  tangential  velocity  at  2000km. 
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Figure  60.  Surface  wave  vertical  velocity  at  2000km. 
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Figure  61.  Surface  wave  radial  velocity  at  2000km,  low-pass  filtered  at  15  sec. 
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Figure  62.  Surface  wave  tangential  velocity  at  2000km,  low-pass  filtered  at  15  sec. 
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Figure  63.  Surface  wave  vertical  velocity  at  2000km,  low-pass  filtered  at  15  sec. 
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4.4.5  Body  Waves 

In  this  section  we  show  far  field  P  and  S  body  waves  from  the  Shoal  and  North  Korea  simulations, 
calculated  for  the  regional  structure  of  Table  5  at  1000  km  from  the  source.  In  each  figure, 
waveforms  are  plotted  for  each  source  depth  along  rows  and  for  azimuths  from  north  counter 
clockwise  to  south  in  steps  of  45  deg.  All  waves  were  computed  for  a  takeoff  angle  of  20°.  P- waves 
for  the  North  Korea  simulations  were  time  advanced  to  approximate  alignment  with  the  main  Shoal 
arrival,  by  a  moveout  equivalent  to  the  travel  time  though  the  NK-Shoal  difference  in  source  depth 
relative  to  the  1-D  surface  level  at  a  slowness  of  l/2000m/s. 
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Figure  64.  Far-field  P-wave  at  1000km  for  Shoal  and  North  Korea  simulations 
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Figure  65.  Far-field  SH-wave  at  1000km  for  Shoal  and  North  Korea  simulations 
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Figure  66.  Far-field  SV-wave  at  1000km  for  Shoal  and  North  Korea  simulations 
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4.5  Decomposition  of  the  Seismic  Source 

Seismic  waves  from  explosions  are  affected  by 
material  properties  at  and  near  the  source,  depth 
and  overburden  pressure,  tectonic  stresses  and 
topography.  We  discussed  topography  in  the 
previous  chapters.  Here  we  look  at  the  effects  of 
depth,  material  properties  and  tectonic  stress 
state.  To  do  this,  we  have  performed  nonlinear 
axisymmetric  calculations  at  13  depths  ranging 
from  150  meters  to  1000  meters,  in  four  earth 
structures  with  distinctly  different  material 
properties,  and  we  have  performed  each 
calculation  with  no  tectonic  stresses  and  with 
tensile  and  compressive  stresses  (Figure  67),  a 
total  of  156  calculations.  Calculations  were 
performed  for  explosion  yields  of  either  10 
kilotons  or  12.5  kilotons.  Normal  containment 
depth  for  these  explosions  is  about  300  meters, 
so  the  calculations  range  from  approximately  Vi 
normal  containment  depth  to  more  than  3  times 
normal  containment  depth.  We  also  performed 
one  dimensional,  spherically  symmetric 
calculations  corresponding  to  each 
axisymmetric  calculation,  using  the  material 
properties  and  overburden  pressure 

corresponding  to  each  depth. 

The  earth  structures  and  material  properties 
used  were  chosen  to  illustrate  the  physical 
phenomena  that  contribute  to  variations  in  near 
source  deformation  and  ground  motion.  These  properties  include  elastic  moduli,  density,  porosity 
and  shear  strength.  Shear  strength  increases  as  a  function  of  pressure  (Figure  68).  Porous  media 
has  a  crush  curve  that  describes  the  change  in  porosity  and  elastic  moduli  as  pores  are  crushed  by 
the  pressure  from  the  explosion.  Of  particular  importance  is  the  reduction  in  strength  after  failure. 
This  is  modeled  in  two  ways:  1)  using  a  damage  model  that  decreases  strength  as  a  function  of 
nonlinear  shear  strain;  and  2)  using  the  effective  stress  model  which  reduces  strength  as  pores  are 
crushed  out.  In  the  effective  stress  model  the  strength  is  determined  by  the  difference  between  the 
pressure  and  the  pore  pressure.  These  become  equal  when  the  pores  are  crushed  out  so  the  strength 
of  the  material  is  determined  by  the  zero  pressure  point  in  Figure  68. 

The  first  two  structures  are  a  Degelen  Mountain  structure  and  a  Climax  Stock  structure.  In  both 
cases  the  material  used  in  the  nonlinear  calculations  was  granite,  with  a  compressional  velocity  of 
5175  m/s,  a  shear  velocity  of  3026  m/s  and  density  of  2600  kg/m3.  The  granite  models  initially 
have  the  same  strength,  however  after  failure  the  Climax  Stock  granite  is  modeled  as  crushed  rock 
with  a  dynamic  coefficient  of  friction  of  0.017,  while  the  Degelen  granite  has  a  post-failure 
coefficient  of  friction  of  0.2,  so  the  Climax  Stock  granite  is  weaker  after  failure.  Strength  reduction 
starts  when  a  cell  reaches  a  nonlinear  strain  of  0.001  and  is  fully  reduced  at  a  nonlinear  strain  of 


Explosion  Locations 


0 


200 

4oo 

8 

o 

<d  600 

Q 

o 

C) 

800 

o 

o 

1000 

o 

Tension  —>• 


•$—  Compression 


200  400 

R(m) 


600 


Figure  67.  Explosions  were  calculated  at  the 
depths  shown  from  150m  to  1000m 
(axisymmetry).  Tectonic  stresses  were 
simulated  by  adding  a  radial  tensile  or 
compressive  stress. 
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0.06.  Shear  modulus  is  also  reduced  after  failure  from  23.8  GPa  to  10  Gpa  for  Climax  Stock  granite 
and  5  GPa  for  Degelen  granite.  Both  of  these  models  were  developed  based  on  modeling  of  near¬ 
field  data  by  Stevens  et  al.  (2003).  In  that  study,  we  had  an  extensive  set  of  near  field  data  from 
Degelen  Mountain  explosions  as  well  as  material  properties  for  Degelen  Mountain  provided  by 
scientists  at  the  Institute  for  the  Dynamics  of  the  Geospheres  in  Moscow,  Russia.  The  Climax 
Stock  model  was  more  recently  used  to  model  the  explosion  Shoal  (Stevens  and  Thompson,  2015), 
so  these  two  models  are  labeled  “Degelen”  and  “Shoal”.  These  models  assume  that  the  pore  space 
is  negligible  and  so  do  not  include  pore  crushing. 


Figure  68.  Strength  of  materials  used  in  this  study.  Right:  expanded  view  of  low  pressure  part 
of  strength  curves.  Strength  is  determined  by  stress  invariants  as  defined  by  Peyton  (1983). 
Strength  shown  in  these  figures  is  approximately  equal  to  shear  stress  divided  by  2.35. 


Specific  Volume  (cm3/gm) 


Figure  69.  Left:  Pahute  Mesa  structure  used  in  depth  of  burial  study  (Day  et  al,  1986)  and  used 
over  a  wider  range  of  depth  in  the  current  study.  Right:  Load/Unload  curves  for  the  three  layers. 
The  other  two  structures  are  based  on  the  depth  of  burial  study  performed  by  Day  et  al.  (1986). 
These  are  models  for  a  three-layered  Pahute  Mesa  structure  and  Shagan  River  granite.  Both  of 
these  models  include  pore  crushing  and  used  the  effective  stress  model  for  strength  reduction.  The 
Pahute  Mesa  structure  and  crush  curves  for  each  layer  are  shown  in  Figure  69.  The  Shagan  River 
structure  was  uniform  granite  with  a  compressional  velocity  of  5018  m/s,  shear  velocity  of  2788 
m/s,  porous  density  of 2700  kg/m3,  and  1%  porosity. 
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4.5.1  Earth  model  with  prestress 

Calculations  were  performed  for  regions  in  both  tensile  and  compressive  prestress.  We  assume 
that  the  water  table  is  close  to  the  surface  and  calculate  the  stress  state  corresponding  to  equilibrium 
for  sliding  on  normal  and  thrust  faults,  respectively  in  a  medium  with  a  coefficient  of  friction  of 
0.6  (Zoback,  2010).  Note  that  this  affects  only  the  horizontal  stress  field;  the  vertical  stress 
corresponds  to  the  weight  of  the  material  above.  The  frictional  equilibrium  state  is  given  by: 


Sl~Pf 

S3-Pf 


P 


+  1  +  //j 


Where  Si  and  S3  are  the  horizontal  and  vertical  principal  stresses  as  discussed  below,  Pf  is  the  pore 
pressure  and  /./  is  the  coefficient  of  friction.  For  q=0.6,  the  expression  on  the  right  hand  side  is 
3.12,  so  the  stress  difference  can  be  quite  large.  In  a  tensile  stress  state  (normal  faulting).  Si  is  Sv, 
the  vertical  stress  and  S3  is  Shmin,  the  minimum  horizontal  stress.  In  a  compressive  stress  state 
(thrust  faulting),  Si  is  Snmax ,  the  maximum  horizontal  stress  and  S3  is  Sv.  Figure  70  shows  the 
calculated  stress  field  as  a  function  of  depth  for  the  granite  models.  Since  we  are  doing 
axisymmetric  calculations,  the  stress  is  uniform  in  all  horizontal  directions  (Sxx=Syy  ±  Szz),  but  the 
horizontal  stress  is  less  than  the  vertical  stress  in  regions  of  extension  (normal  faulting)  and  greater 
than  the  vertical  stress  in  regions  of  compression  (thrust  faulting).  Szz  is  equal  to  Sv,  the  vertical 
stress  determined  by  the  weight  of  the  material  above. 


Stress  (MPa)  Stress  (MPa) 

Figure  70.  Left:  tensile  stress  state  corresponding  to  a  state  in  balance  with  normal  faulting  with 
a  coefficient  of  friction  of  0.6.  Right:  compressive  stress  state  in  balance  with  thrust  faulting  with 
a  coefficient  of  friction  of  0.6.  Dashed  lines  are  the  limiting  stress;  solid  lines  are  the  stress  used 
in  the  calculations. 


4.5.2  Near  field  nonlinear  deformation 

Shallow  explosions  generate  a  nearly  conical  region  of  nonlinear  deformation  around  the  explosion 
and  up  to  the  surface  while  deeper  explosions  generate  a  more  spherical  region  of  nonlinear 
deformation  around  the  explosion  plus  some  near-surface  spall.  Figure  71  through  Figure  74  show 
the  region  of  nonlinear  deformation  for  the  four  structures  at  depths  of  150,  300,  600  and  1000 
meters. 


Approved  for  public  release;  distribution  is  unlimited. 

67 


-400 


Yield  -  at  2.00s 


Yield  -  at  2.00s 


^300 

-200 


300' 


400  — 

0  200  400  600 

Radial  Distance  (m) 


-400 

^300 


300 

400  — 

0  200  400  600 

Radial  Distance  (m) 


Figure  71.  Region  of  nonlinear  deformation  for  an  explosion  at  150  meters  depth.  Top  left: 
Degelen;  top  right:  Shoal;  bottom  left:  Shagan;  bottom  right:  Pahute.  The  Degelen  explosion 
was  10  kt;  the  other  three  12.5  kt. 
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Figure  72.  Region  of  nonlinear  deformation  for  an  explosion  at  300  meters  depth.  Top  left: 
Degelen;  top  right:  Shoal;  bottom  left:  Shagan;  bottom  right:  Pahute. 
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Figure  73.  Region  of  nonlinear  deformation  for  an  explosion  at  600  meters  depth.  Top  left: 
Degelen;  top  right:  Shoal;  bottom  left:  Shagan;  bottom  right:  Pahute. 
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Figure  74.  Region  of  nonlinear  deformation  for  an  explosion  at  1000  meters  depth.  Top  left: 
Degelen;  top  right:  Shoal;  bottom  left:  Shagan;  bottom  right:  Pahute. 
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4.5.3  Near  field  nonlinear  deformation  in  a  tectonic  stress  field 

Figure  75  to  Figure  77  show  the  near  field  nonlinear  deformation  for  explosions  at  depths  of  150m, 
300m  and  600m  in  compressive  and  tensile  stress  fields.  In  all  cases  there  is  more  deformation  in 
tensile  stress  fields  and  less  in  compressive  stress  fields.  There  are  two  effects  operating  here:  1) 
both  compressive  and  tensile  stress  fields  increases  the  shear  stress;  and  2)  the  compressive  stress 
field  increases  the  pressure  thereby  increasing  strength  while  the  tensile  stress  field  decreases  the 
pressure  and  so  decreases  strength. 
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Figure  75.  Region  of  nonlinear  deformation  for  an  explosion  in  a  prestress  field  at  150  meters 
depth.  Left  figures  are  tensile  stress  fields,  right  figures  are  compressive  stress  fields.  Top 
row:  Degelen;  second  row:  Shoal;  third  row:  Shagan;  bottom  row:  Pahute. 
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Figure  76.  Region  of  nonlinear  deformation  for  an  explosion  in  a  prestress  field  at  300  meters 
depth.  Left  figures  are  tensile  stress  fields,  right  figures  are  compressive  stress  fields.  Top 
row:  Degelen;  second  row:  Shoal;  third  row:  Shagan;  bottom  row:  Pahute. 
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Figure  77.  Region  of  nonlinear  deformation  for  an  explosion  in  a  prestress  field  at  600  meters 
depth.  Left  figures  are  tensile  stress  fields,  right  figures  are  compressive  stress  fields.  Top 
row:  Degelen;  second  row:  Shoal;  third  row:  Shagan;  bottom  row:  Pahute. 
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4.5.4  Propagation  using  the  representation  theorem 

The  representation  theorem  was  used  to  calculate  surface  waves  at  2000  km,  regional  phases  at 
250  km,  and  far-field  body  waves  from  each  calculation.  An  equivalent  point  source  calculation 
was  made  for  each  using  the  reduced  velocity  potential  (RVP)  convolved  with  a  plane-layered 
Green’s  function  for  each  phase.  These  were  compared  to  determine  the  effect  of  the  nonlinear 
interaction  with  the  free  surface  compared  to  an  elastic  calculation.  The  earth  structures  used  for 
calculations  of  regional  phases  and  surface  waves  are  shown  in  Figure  78.  These  are  regional 
structures  for  each  site  with  the  material  properties  of  the  nonlinear  calculations  in  the  upper  part 
of  the  structure.  For  the  Pahute  structure,  the  two  slow,  shallow  layers  were  removed  and  the 
properties  of  the  third  layer  used  all  the  way  to  the  surface.  The  nonlinear  calculations  also  merged 
into  this  structure  outside  of  the  source  region.  The  Degelen  and  Shagan  structures  used  the  earth 
model  for  that  region  derived  from  surface  waves  by  Stevens  et  al  (2008).  They  differ  only  in  the 
material  properties  of  the  top  layer,  which  is  lower  velocity  for  Shagan  than  for  Degelen.  The 
Pahute  and  Shoal  structures  are  derived  from  an  earth  model  for  the  Nevada  Test  Site. 


Vp  (km/s)  Vs  (km/s) 

Figure  78.  Vp  and  Vs  vs.  depth  for  the  four  structures  used  in  the  calculations. 


Approved  for  public  release;  distribution  is  unlimited. 

74 


4.5.5  Surface  wave  results 


Surface  waves  were  calculated  as  discussed  above  for  each  axisymmetric  calculation  in  each  of 
the  four  structures,  and  also  for  the  equivalent  point  source  calculation  at  each  depth.  We  first 
summarize  the  results  in  Figure  79  and  Figure  80  and  then  show  some  waveform  examples.  We 
find  that  surface  waves  from  one-dimensional  spherically  symmetric  calculations  and  two- 
dimensional  axisymmetric  calculations  agree  fairly  well  and  show  a  slow  decrease  with  depth  for 
depths  greater  than  about  300  meters,  however  at  shallower  depths  the  free  surface  interaction 
causes  a  sharp  decrease  in  the  surface  wave  amplitude.  Surface  waves  are  reduced  by  compressive 
tectonic  strain  release  and  increased  by  tensile  tectonic  strain  release,  with  compressive  being  the 
stronger  effect.  The  shallower  calculations  in  the  Climax  stock  (Shoal)  granite  calculation,  with 
compressive  tectonic  strain  release  are  polarity  reversed.  Surface  waves  from  the  Pahute 
calculation  are  more  complicated  because  of  the  high  porosity  unsaturated  media  that  reduces 
amplitudes  for  depths  less  than  580  meters. 
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Figure  79.  20  second  surface  wave  magnitudes  from  all  axisymmetric  and  spherically  symmetric 
calculations.  Left:  Degelen  granite;  Right:  Climax  stock  granite.  Black  circles  indicate  surface 
waves  with  polarity  reversals. 
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Figure  80.  20  second  surface  wave  magnitudes  from  all  axisymmetric  and  spherically  symmetric 
calculations.  Left:  Shagan  River  granite;  Right:  Pahute  Mesa  Rhyolite. 
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Figure  81  and  Figure  82  show  the  surface  waves  for  the  point  source  calculation  and  the 
axisymmetric  calculation  at  explosion  depths  of  150  and  300  meters  in  the  four  structures.  All  are 
reduced  in  amplitude  substantially  at  150  meters,  and  the  Pahute  calculation  while  not  actually 
reversed,  is  clearly  phase  shifted.  At  300  meters  the  point  source  and  axisymmetric  calculations 
are  much  more  similar,  particularly  for  the  Degelen  and  Shagan  models.  Amplitudes  are  slightly 
reduced  by  the  free  surface  interaction  for  the  Shoal  and  Pahute  models. 

A  reduction  in  surface  wave  amplitude  due  to  near-surface  damage  was  proposed  by  Patton  and 
Taylor  (2008)  and  Patton  (2016).  These  calculations  support  that  conclusion  for  shallow  explosion 
sources,  with  the  reduction  starting  at  about  normal  containment  depth.  However  tectonic  strain 
release  also  has  a  strong  effect  that  can  both  increase  and  decrease  surface  wave  amplitudes  for 
tensile  and  compressive  stresses,  respectively.  Increasing  depth  also  decreases  surface  wave 
amplitudes  due  to  increasing  overburden  pressure,  so  the  combination  of  the  near-surface  decrease 
and  depth  decrease  causes  a  peak  surface  wave  amplitude  at  a  depth  near  normal  containment 
depth. 
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Figure  81.  Surface  waves  for  an  explosion  depth  of  150  meters  in  four  structures.  Figures 
show  point  source  seismograms  in  blue  and  axisymmetric  calculation  without  prestress  in 
blue.  Top  left:  Degelen,  top  right:  Shoal,  lower  left:  Shagan,  lower  right:  Pahute. 
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Figure  82.  Surface  waves  for  an  explosion  depth  of  300  meters  in  four  structures.  Figures 
show  point  source  seismograms  in  blue  and  axisymmetric  calculation  without  prestress  in 
blue.  Top  left:  Degelen,  top  right:  Shoal,  lower  left:  Shagan,  lower  right:  Pahute. 


Surface  waves  are  generated  primarily  by  the  static  horizontal  displacements  caused  by  the 
explosion.  Figure  83  -  Figure  85  show  the  final  horizontal  displacement  at  depths  of  200,  370  and 
800  meters,  respectively,  in  the  Degelen  and  Shoal  structures  and  in  tensile,  no  and  compressive 
prestress.  In  all  cases  tensile  stress  increases  the  horizontal  displacement  and  compressive  stress 
decreases  it.  In  the  200  meter  Shoal  case  with  compressive  stress,  the  horizontal  displacement  is 
negative,  leading  to  the  Rayleigh  wave  reversals  shown  in  Figure  79. 
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Figure  83.  Final  horizontal  displacement  at  200m  depth.  Left:  tensile  prestress.  Center:  no 
prestress.  Right:  compressive  prestress.  Top:  Degelen;  bottom:  Shoal.  The  horizontal 
displacement  with  Shoal  in  compression  is  almost  entirely  negative,  leading  to  the  Rayleigh 
wave  reversals  shown  in  Figure  79. 
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Figure  84.  Final  horizontal  displacement  at  370m  depth.  Left:  tensile  prestress.  Center:  no 
prestress.  Right:  compressive  prestress.  Top:  Degelen;  bottom:  Shoal. 
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Figure  85.  Final  horizontal  displacement  at  800m  depth.  Left:  tensile  prestress.  Center:  no 
prestress.  Right:  compressive  prestress.  Top:  Degelen;  bottom:  Shoal. 

To  find  out  why  the  near  surface  amplitude  reduction  occurs,  we  ran  an  additional  set  of 
calculations  in  the  Degelen  model,  but  with  the  strength  increased  enough  that  no  nonlinear 
deformation  occurred  (referred  to  as  the  “Hard”  calculation),  so  the  problem  is  purely  elastic.  We 
then  compared  the  results  of  that  calculation  with  the  Degelen  and  Shoal  calculations. 
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Figure  86.  Left:  Ms  from  ID  and  2D  “Hard”  calculation.  Middle:  Hard  model  radial 
displacement.  Right:  Degelen  Model  radial  displacement.  Both  for  150  meters  depth. 

Figure  86  shows  the  calculated  surface  wave  amplitudes  for  the  hard  model  and  the  radial 
displacement  for  the  Hard  and  Degelen  models,  with  the  scale  on  the  hard  model  increased  by  a 
factor  of  22.  Now  as  expected  the  surface  waves  show  no  depth  dependence  for  either  the  ID  or 
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2D  models.  The  patterns  of  radial  displacement  are  very  similar,  but  there  is  a  reduction  in 
horizontal  displacement  that  increases  toward  the  surface  in  the  Degelen  model.  We  can  see  this 
more  clearly  by  examining  the  difference  between  the  two  final  displacement  fields,  which  are 
shown  in  Figure  87.  Above  the  source  the  displacement  difference  is  upward  and  outward  where 
the  material  is  pushed  up  and  away  from  the  explosion.  However  a  short  distance  out  the 
displacement  difference  changes  to  inward  and  downward,  almost  as  if  there  were  a  thrust  fault 
extending  from  the  explosion  to  the  surface.  This  is  similar  to  the  mechanism  proposed  by  Masse 
(1981),  who  suggested  that  induced  thrust  faulting  above  the  explosion  was  responsible  for 
Rayleigh  Wave  reversals  and  other  effects  attributed  to  tectonic  strain  release. 
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Figure  87.  Radial  (left)  and  Vertical  (right)  displacement  field  difference  between  the  Degelen 
and  (scaled)  Hard  models  for  an  explosion  at  150  meters  depth. 
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Figure  88.  Vectors  of  displacement  difference.  Left:  Degelen  at  150m  depth.  Right:  Shoal  at 
250m  depth. 
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The  shallowest  explosions  appear  to  be  quite  consistent  with  the  Masse  model,  but  deeper 
explosions  have  more  complex  displacement  fields.  Figure  88  shows  vectors  of  displacement 
difference  between  the  Degelen  and  Hard  explosions  at  150  meters  depth  and  between  the  Shoal 
and  Hard  explosions  (scaled  by  27.5).  At  250  meters  there  is  a  “whirl”  of  displacement  to  the  right 
of  the  cavity  and  inward  motion  along  the  surface. 

4.5.6  Full  Wave  Regional  Waveforms 

We  calculated  full  waveform  seismograms  for  all  of  the  calculations  at  a  distance  of  250  km  using 
the  earth  models  described  earlier.  Example  seismograms  from  the  source  at  300  meters  depth 
(without  tectonic  release)  are  shown  in  Figure  89.  Each  figure  shows  a  waveform  from  the 
nonlinear  calculations  propagated  using  the  representation  theorem,  and  also  a  waveform  from  a 
point  source  plane-layered  medium  convolved  with  the  RVP  from  the  corresponding  spherically 
symmetric  calculation. 
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Figure  89.  Example  full  waveform  seismograms  for  the  four  earth  models.  Seismograms  labeled 
2D  Axisymmetric  are  generated  from  the  nonlinear  axisymmetric  calculations  using  the 
representation  theorem.  Seismograms  labeled  ID  RVP  are  point  source  plane-layered  full 
waveform  seismograms  convolved  with  an  RVP  for  the  source. 

We  are  interested  in  how  the  nonlinear  interaction  with  the  free  surface  and  tectonic  release  affect 
regional  phase  amplitudes.  To  do  this,  we  filtered  each  seismogram  with  a  3  pole  two-pass 
Butterworth  filter  with  a  passband  of  0.6  to  3.0  Hz,  then  found  the  peak  amplitude  within  arrival 
time  windows  of  35-55  seconds  (4.5-7. 1  km/s)  for  P-waves,  55-65  seconds  (3. 8-4.5  km/s)  for  S- 
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waves,  and  65-78  seconds  (3.2-3. 8  km/s)  for  Lg.  Measurements  for  all  calculations  are  shown  in 
the  following  figures. 

Figure  90  shows  regional  P-wave  amplitudes  for  all  calculations.  There  is  a  slow  decrease  in 
amplitude  with  depth  caused  by  the  increasing  overburden  pressure.  The  2D  and  ID  calculations 
for  the  three  granite  structures  are  very  similar  and  tectonic  stresses  appear  to  have  little  effect  on 
the  P-waves.  The  Pahute  structure  also  shows  little  effect  on  P-waves  from  tectonic  stresses.  The 
ID  Pahute  calculations  have  higher  amplitudes  than  the  2D  calculations,  however  the  ID 
calculations  were  embedded  in  the  structure  without  the  two  surface  layers,  so  the  ID  and  2D 
Pahute  calculations  may  not  be  directly  comparable. 


Figure  90.  P-wave  amplitudes  for  each  structure.  Blue  points  are  the  ID  calculation,  Red  points 
are  the  2D  axisymmetric  calculations  with  and  without  tectonic  release. 

Figure  91  shows  regional  S-wave  amplitudes  for  all  calculations.  Now  the  effect  of  the  nonlinear 
interaction  with  the  free  surface  is  much  more  apparent.  In  all  of  the  granite  structures  at  shallower 
depths  the  surface  interaction  generates  substantially  more  S-waves  compared  to  the  point  source, 
which  only  generates  S-waves  through  pS  reflection. 
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Figure  91.  S-wave  amplitudes  for  each  structure.  Blue  points  are  the  ID  calculation,  Red  points 
are  the  2D  axisymmetric  calculations  with  and  without  tectonic  release. 

Figure  92  shows  Lg  amplitudes  for  all  calculations  (with  the  caveat  that  any  enhancement  of 
Lg  due  to  scattering  is  not  included).  The  free  surface  interaction  has  surprisingly  little  effect  on 
Lg,  with  the  Degelen  and  Shoal  structures  showing  a  small  near-surface  decrease  similar  to 
the  fundamental  mode  surface  waves.  For  deeper  sources,  tectonic  release  does  affect  the 
Lg  amplitudes,  and  does  so  opposite  the  effect  on  fundamental  mode  surface  waves  with 
compression  increasing  Lg  and  tension  decreasing  it. 
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Figure  92.  Lg  amplitudes  for  each  structure.  Blue  points  are  the  ID  calculation,  Red  points  are 
the  2D  axisymmetric  calculations  with  and  without  tectonic  release. 
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4.5.7  Far-field  Body  Waves 

Far-field  P  and  SV  waves  were  calculated  for  all  of  the  nonlinear  calculations  using  the 
representation  theorem.  Because  we  are  interested  in  the  effects  of  the  nonlinear  interaction  with 
the  free  surface  rather  than  crustal  reflections,  we  use  only  the  part  of  each  structure  above  1530 
meters  and  calculate  the  P  and  SV  waves  as  they  leave  the  base  of  the  structure.  P  waves  are 
calculated  at  a  takeoff  angle  of  10°,  SV  waves  at  20°.  Example  waveforms  for  a  depth  of  250  meters 
are  shown  in  Figure  93.  The  main  physical  effect  operating  here  is  the  degradation  of  the  pP  and 
pS  phases  by  the  nonlinear  interaction  with  the  free  surface.  The  initial  P-wave  amplitude  is  very 
similar  between  the  ID  and  2D  calculations  because  the  initial  downgoing  P-wave  is  unaffected 
by  the  free  surface,  however  subsequent  reflected  phases  are  decreased  in  amplitude,  distorted  and 
spread  out  in  time.  The  Pahute  SV  waveforms  are  complex  because  of  reflections  within  the  low 
velocity  surface  layers.  The  waveforms  shown  in  Figure  93  are  low-pass  filtered  at  10  Hz,  and  in 
the  frequency  band  the  pP  tends  to  increase  the  P-wave  amplitude.  That  combined  with  the  larger 
RVP  for  shallow  events  causes  the  ID  amplitudes  to  increase  at  shallow  depths,  while  the  P-wave 
amplitudes  from  the  2D  calculations  remain  almost  constant  with  depth. 


0  0.5 


1.5  2 


0  0.5 


1.5  2 


Time  (s)  Time  (s) 


0  0.5 


1.5  2 


0  1  2  3  4  5 


Time  (s)  Time  (s) 

Figure  93.  Far-field  P  and  SV  waves  from  ID  and  2D  calculations  at  a  depth  of  250  meters  in  the 
four  earth  structures.  No  tectonic  stresses.  Mark  shows  measured  amplitude  for  the  2D 
waveform.  All  waveforms  are  low-pass  filtered  at  10  Hz. 


Figure  94  and  Figure  95  show  the  amplitude  measurements  for  all  P  and  SV  waveforms  for  all 
depths,  with  and  without  tectonic  stresses  in  the  four  structures. 
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Figure  94.  Far  field  P-wave  amplitudes  for  each  structure.  Blue  points  are  the  ID  calculation, 
Red  points  are  the  2D  axisymmetric  calculations  with  and  without  tectonic  release. 


Figure  95.  Far  field  SV-wave  amplitudes  for  each  structure.  Blue  points  are  the  ID 
calculation,  Red  points  are  the  2D  axisymmetric  calculations  with  and  without  tectonic 
release. 
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Figure  95  (continued).  Far  field  SV-wave  amplitudes  for  each  structure.  Blue  points  are  the  ID 
calculation,  Red  points  are  the  2D  axisymmetric  calculations  with  and  without  tectonic  release. 
The  results  described  above  change  at  lower  frequencies  as  the  direct  P  and  pP  destructively 
interfere  at  low  frequencies.  The  effect  of  this  is  to  reduce  the  P-wave  amplitudes  in  the  ID 
calculations  that  increased  at  the  higher  frequencies,  making  the  ID  and  2D  amplitudes  more 
consistent  and  almost  independent  of  depth.  Example  waveforms  are  shown  in  Figure  96. 
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Figure  96.  Same  waveforms  as  Figure  93  lowpass  filtered  at  1  Hz. 

Figure  97  and  Figure  98  show  the  amplitude  measurements  for  all  P  and  SV  waveforms  filtered 
at  1  Hz  for  all  depths,  with  and  without  tectonic  stresses  in  the  four  structures 
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Figure  97.  Far  field  P-wave  amplitudes,  low-pass  filtered  at  1  Hz,  for  each  structure.  Blue  points 
are  the  ID  calculation,  Red  points  are  the  2D  axisymmetric  calculations  with  and  without 
tectonic  release. 


Figure  98.  Far  field  SV-wave  amplitudes,  low-pass  filtered  at  1  Hz,  for  each  structure.  Blue 
points  are  the  ID  calculation,  Red  points  are  the  2D  axisymmetric  calculations  with  and 
without  tectonic  release. 
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Figure  98  (continued).  Far  field  SV-wave  amplitudes,  low-pass  filtered  at  1  Hz,  for  each 
structure.  Blue  points  are  the  ID  calculation,  Red  points  are  the  2D  axisymmetric  calculations 
with  and  without  tectonic  release. 


4.5.8  Conclusions  from  Seismic  Source  Decomposition  Analysis 

Figure  99  shows  the  Ms  measurements  from  East  Kazakh  explosions  and  calculations  of  Ms  with 
and  without  compressive  and  tensile  tectonic  release.  The  East  Kazakh  Ms  values  are  in  good 
agreement  with  the  calculations,  particularly  with  compressive  tectonic  release.  The  North  Korean 
Ms  values  are  still  very  high  relative  to  the  calculations,  however,  even  for  tensile  tectonic  release. 
Assuming  the  yields  are  correct  for  these  events  (they  are  estimates  derived  from  P-waves),  there 
are  several  possible  explanations  for  this.  As  noted  in  the  previous  sections  on  3D  calculations 
with  topography,  there  is  an  amplification  effect  for  explosions  near  the  base  of  a  mountain,  and 
the  North  Korean  events  all  appear  to  be  near  the  base  of  a  mountain.  Another  possible  explanation 
is  underestimation  of  the  effect  of  tectonic  strain  release.  The  tectonic  stresses  were  derived  from 
equilibrium  with  frictional  sliding,  however  the  granite  strength  model  is  stronger  than  frictional 
sliding,  so  accounting  for  frictional  slip  would  increase  the  effect  of  tectonic  release.  Finally,  as 
noted  in  Figure  100,  a  weaker  granite  model  such  as  that  used  for  modeling  Shoal,  generates  larger 
surface  waves  than  the  stronger  Degelen  model. 
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Figure  99.  Ms-log(Yield)  vs.  scaled  depth.  Left:  Measurements  from  East  Kazakh  and  North 
Korean  Data.  Right:  Calculations  for  Degelen  structure. 
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Figure  100.  Calculations  for  Shoal  structure.  Black  circles  indicate  polarity-reversed  Rayleigh 
waves. 

5.  Improvements  to  CRAM3D  Codes 
Accelerated  Equilibrium  Calculation 

The  longest  part  of  a  CRAM3D  simulation  with  topography  is  the  computation  of  an  equilibrium 
starting  grid.  There  are  two  reasons  for  this.  First,  uncompensated  topographic  effects  that  are  left 
for  the  simulation  to  resolve  are  large  in  space  and  mass,  resulting  in  correspondingly  large 
oscillatory  grid  motions  that  take  proportionally  longer  to  die  out,  up  to  several  times  more 
simulation  time  than  for  most  of  the  energy  from  an  explosion  to  exit  or  be  absorbed  by  the  grid. 
Second,  numerical  stability  requires  the  size  of  simulation  time  steps  to  be  proportional  to  the 
smallest  separation  between  any  two  nodes  in  the  grid.  The  nodes  in  the  inner  grid  at  the  source 
cavity  boundary  are  typically  ~0.5  m  apart,  compared  to  O(10m)  in  the  outer  grid.  In  the  first  l/10th 
second  of  the  explosion  simulation,  the  cavity  can  grow  by  many  times,  separating  those  nodes 
laterally,  and  after  the  main  shock  has  passed  we  can  merge  layers  of  zones  around  the  cavity, 
thereby  increasing  the  minimum  node  separation  and  speeding  the  rest  of  the  simulation.  Neither 
applies  to  establishing  the  equilibrium,  where  the  cavity  and  grid  must  remain  intact. 

To  speed  the  equilibrium  computations,  we  have  implemented  the  ability  to  compute  and  apply 
topographic  overburden  pressures  that  are  in  isostatic  equilibrium.  This  is  based  on  an 
axisymmetric  solution  for  gravitational  equilibrium  beneath  a  sloping  surface  (Stevens  et  al,  1990) 
that  also  provides  good  equilibrium  in  3D.  In  this  solution,  the  topography  above  the  base  surface 
level  (that  of  the  plane-layer  earth  surrounding  the  grid)  contributes  only  1/3  to  the  overburden 
pressure.  The  contributions  are  easily  summed  for  the  outer  grid.  The  inner  grid  state  is 
approximated  by  that  of  a  full  outer  grid.  As  discussed  earlier,  we  also  break  the  equilibrium 
calculation  into  two  parts,  first  with  the  inner  grid  replaced  with  outer  grid  zones  so  that  we  do  not 
have  the  tiny  zone  problem  discussed  above,  followed  by  a  shorter  calculation  with  the  inner  grid 
included. 

In  addition,  we  have  improved  the  process  by  which  motions  are  damped  out  during  the 
equilibrium  simulation.  In  the  past,  we  have  reduced  grid  velocity  by  a  small  factor  at  each  time 
step,  with  the  goal  of  short  circuiting  large  grid  oscillations  that  would  otherwise  result  if  starting 
with  a  large  uncompensated  topographic  mass.  But  the  schemes  for  selecting  the  damping  have 
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been  rudimentary  and  their  performance  not  well  characterized.  We  have  adopted  a  new  scheme 
wherein  the  net  grid  motion  is  calculated  at  each  step  and  damping  only  applied  when  the  grid  is 
already  slowing,  i.e.,  the  net  acceleration  is  opposite  the  net  velocity.  We  performed  simulations 
with  large  uncompensated  topographic  loads  and  varied  the  damping  factor  to  characterize  the 
scheme.  The  results  are  illustrated  in  Figure  101.  Each  plots  the  maximum  grid  velocity  with 
simulation  time  for  a  given  value  of  the  damping.  The  green  curve  is  undamped,  displaying  grid 
oscillations  that  continue  for  many  simulation  seconds.  In  contrast,  grid  motion  is  down  to  an 
almost  negligible  level  after  only  2  seconds  of  simulation  time  with  the  additional  damping. 

With  these  improvements,  we  have  reduced  equilibrium  simulation  times  by  a  factor  of  as  much 
as  5.  For  runs  that  previously  took  several  days  to  complete,  this  is  a  very  significant  change. 


Time  to  equilibrium  with  large  topography 


Figure  101.  The  effects  of  a  dynamically  varied  damping  in  reducing  simulation  time  for 
attaining  grid  equilibrium.  The  green  curve  is  undamped. 
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6.  Conclusions 


We  retrieved  data  from  small  Semipalatinsk  (Shagan  River  and  Degelen)  explosions  at  Borovoye 
and  measured  Ms  for  events  as  small  as  1.4  kilotons.  We  find  that  Ms:yield  is  very  consistent  with 
Ms  measured  from  larger  events,  with  a  slope  of  1.0,  and  very  consistent  with  the  global  Ms:yield 
distribution  curve  (Stevens  and  Murphy,  2001),  and  they  are  very  inconsistent  with  Ms  vs. 
estimated  yield  for  the  three  North  Korean  explosions,  which  lie  almost  a  full  magnitude  unit  above 
the  Semipalatinsk  curve.  This  shows  that  the  high  Ms  of  the  North  Korean  explosions  is  not  the 
result  of  a  nonunity  Ms:yield  curve,  and  that  there  is  something  very  different  about  generation  of 
surface  waves  from  the  North  Korean  explosions  and  the  Semipalatinsk  explosions.  We  have 
previously  hypothesized  that  this  difference  could  be  due  to  differences  tectonic  stress  state  and 
tectonic  release  by  the  explosions  (Murphy  et  al,  2013),  and  this  remains  a  likely  explanation.  We 
also  find  in  this  study  an  amplification  of  surface  waves  for  explosions  near  the  base  of  a  mountain 
with  steep  topography. 

We  examined  the  effect  of  depth  beneath  a  mountainside  or  other  steep  slope  on  seismic  waves  by 
performing  three  3D  nonlinear  calculations  of  an  explosion  beneath  a  steep  slope.  The  first,  Slope 
1,  is  placed  beneath  the  mountain;  the  second,  Slope  2,  is  higher  and  inside  the  mountain;  while 
the  third,  Slope  3,  is  inside  the  mountain  with  the  yield  increased  by  a  factor  of  4.  We  then  did  four 
additional  3D  calculations  using  the  actual  topography  of  the  North  Korean  test  site. 

We  find  that  generation  of  long  period  surface  waves  is  strongly  dependent  on  location  relative  to 
the  slope.  Long  period  surface  waves  are  generated  almost  entirely  by  the  static  horizontal 
displacement  in  the  upper  few  hundred  meters  below  the  surface.  If  the  explosion  is  inside  the 
mountain  and  close  enough  to  the  free  surface  that  stress  is  relieved  by  the  sides  of  the  mountain, 
then  subsurface  horizontal  displacement  is  reduced  and  surface  wave  amplitudes  are  reduced 
dramatically.  However,  for  explosions  that  are  near  the  base  of  the  mountain,  long  period  surface 
waves  are  increased  in  amplitude  by  as  much  as  a  factor  of  2. 

In  contrast,  regional  P  and  S  phases  are  slightly  larger  for  the  explosion  in  the  mountain  than 
beneath  it,  especially  for  paths  that  go  across  the  mountain.  Otherwise,  there  is  little  difference  in 
regional  phases  between  the  flat  (Shoal)  and  topographic  (Slope  1  and  Slope  2)  calculations.  The 
North  Korean  calculations  show  variability  in  regional  P  and  S  phase  amplitudes  as  a  function  of 
azimuth  and  as  a  function  of  depth,  however  the  average  changes  are  small. 

For  the  Slope  calculations,  far  field  P  body  waves  are  slightly  enhanced  by  the  presence  of  the 
mountain  due  to  reduction  of  the  pP  interference.  That  is,  the  pP  phase  with  the  mountain  present 
interferes  less  with  the  direct  P  than  it  does  for  a  flat  surface.  For  Slope  1  this  occurs  primarily  for 
the  cross-slope  direction,  while  for  Slope  2  it  occurs  in  all  directions.  Far  field  SV  body  waves  are 
reduced  in  amplitude  relative  to  the  flat  case  because  the  pS  phase  which  is  very  strong  for  a  flat 
surface  is  reduced  in  amplitude  for  the  mountain.  Far  field  SH  waves  which  do  not  exist  for  the 
flat  surface  are  comparable  in  amplitude  to  SV  for  the  mountain  cases.  Results  for  the  North  Korea 
calculations  are  similar  except  that  the  P-wave  amplitudes  are  more  consistent  and  show  minimal 
effect  from  either  topography  or  source  depth.  SV  waves  are  also  more  consistent  than  in  the  Slope 
calculations,  probably  because  the  NK  mountain  is  larger  than  the  mountain  used  in  the  Slope 
calculations  and  the  surface  variation  occurs  over  a  larger  area. 
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In  addition  to  the  3D  calculations,  we  performed  a  large  number  of  2D  axisymmetric  calculations 
for  four  earth  structures,  with  and  without  tectonic  release  and  covering  a  range  of  depths  from 
150m  to  1000m  to  examine  the  effects  of  depth,  material  properties  and  tectonic  stress  state.  We 
find  that  P-waves  are  relatively  insensitive  to  these  effects  because  the  initial  downgoing  P-wave 
is  unaffected  by  interaction  with  the  free  surface.  Surface  waves,  however,  are  strongly  affected 
by  all  of  these  effects.  Surface  waves  decrease  in  amplitude  at  shallow  depths  where  there  is  strong 
interaction  with  the  free  surface,  and  they  are  decreased  by  compressive  tectonic  release.  Surface 
waves  are  also  decreased  by  increasing  depth.  Tensile  tectonic  stresses  amplify  surface  waves, 
although  the  effect  of  compressive  stresses  to  decrease  surface  waves  is  stronger.  Surface  waves 
are  largest  for  explosions  at  or  slightly  below  normal  containment  depth  where  the  interaction  with 
the  free  surface  is  moderate  and  the  increasing  pressure  with  depth  does  not  reduce  them.  This 
effect  is  amplified  by  topography  where  surface  waves  from  explosions  within  the  mountain  are 
reduced  for  the  reasons  explained  earlier,  and  explosions  near  the  base  of  the  mountain  are 
increased  because  the  pressure  and  horizontal  stress  are  less  than  at  comparable  depth  in  a  region 
without  topography.  The  North  Korean  events  appear  to  be  optimally  located  for  surface  waves: 
they  are  near  the  base  of  the  mountain  and  in  a  region  believed  to  have  tensile  tectonic  stress. 
However,  even  under  these  optimum  conditions,  calculations  do  not  produce  surface  waves  quite 
as  large  as  those  observed.  A  possible  explanation  is  underestimation  of  the  effect  of  tectonic  strain 
release.  The  tectonic  stresses  were  derived  from  equilibrium  with  frictional  sliding,  however  the 
granite  strength  model  is  stronger  than  frictional  sliding,  so  accounting  for  frictional  slip  would 
increase  the  effect  of  tectonic  release. 
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Appendix  -  CRAM3D  Update  and  Deliverable 


These  are  Release  Notes  for  the  Leidos  2017-05  release  of  source 
code  for  the  CRAM3D  and  related  programs  under  AFRL  contract 
FA9453-14-C-0257 

CHANGES : 

CRAM3D 


1 .  Added  Parameters : 
iso_comp  (switch) 

A  factor  between  0  and  1  that  specifies  the  degree  to 
which  the  initial  stresses  account  for  the  topography 
isostatically  vs.  material  strength 

2.  Updated  initial  pressure  caclulation  with  isostatic 
topographic  loading 

-  1/3  weighted  topo  load  in  mu  recursion 

-  axial  stresses  in  outer/inner  grid 

-  stresses  in  inner  grid  are  not  accurate  with  multiple 
materials 

-  also  allows  for  dynamic  compensation  also 

3.  Improved  damping  scheme  for  equilibrium  run 

-  application  of  damping  based  on  v_avg* (dv/dt ) _avg  over  grid 

4.  Improved  process  logging. 

5.  Saves  zone  masses  to  disk 

6 .  cram_plot_dump 

-  adds  ability  to  plot  vertical  stress,  radial  displacement, 
volumetric  strain 

-  simplified  usage 

-  fix  for  new  colorbar  behavior  in  R2015a 

-  replaces  pcolor  with  surf  in  response  to  new  behavior  in 
R2015a 

-  fixes  symbol  plotting  by  replacing  scatter  with  scatter3 

-  Improved  exception  handling  in  various  of  plotting  routines. 


INSTALLATION: 

Untar  dist.tar  in  an  empty  target  directory  and  follow  the 
build/test  instructions  in  the  README 
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AFRL 

AFSPC 

AFWA 

NPE 

TOA 


List  of  Symbols,  Abbreviations,  and  Acronyms 

Air  Force  Research  Laboratory 
Air  Force  Space  Command 
Air  Force  Weather  Agency 
Non-proliferation  Experiment 
Take-off  Angle 
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